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Summary 



In this thesis, I study hydrodynamical models of shm accretion disks — advective, 
optically thick disks which generalize the standard models of radiatively efficient thin 
disks to all accretion rates. I start with a general introduction to the theory of accretion 
onto compact objects. It is followed by a derivation of the commonly-used standard 
models of thin disks. In the subsequent section I introduce the equations describing 
slim disks, explain the numerical methods I used to solve them and discuss properties of 
such solutions. I also give a general derivation of non-stationary equations and present 
the time evolution of thermally unstable accretion disks. I introduce a state-of-the-art 
approach coupling the radial and vertical structures of an advective accretion disk and 
discuss the improvements it brings to vertically-averaged solutions. I also present a 
numerical model of self-illuminated slim accretion disks. Finally, I present and discuss 
applications of slim accretion disks: estimating of spin of the central black hole in LMC 
X-3 through X-ray continuum fitting basing on high-luminosity data, spinning-up of 
black holes by super-critical accretion fiows and normalizing of magnetohydrodynamical 
global simulations. 



This thesis presents and extends results included in the following papers: 

S4DOWSKI (2009) - Sections 4.I — 4.3, 
Abramowicz et AL. (2010) — Section 4-4, 
XUE ET AL. (2011) - Chapter 5, 
S4DOWSKI ET AL. (2011) - Chapter 6, 
S4DOWSKI ET AL. (2011) — Section 8.1. 

Unpublished material is included in Section 4-5.2 and in Chapter 7. 



Chapter 1 
Introduction 



Microquasars and Active Galactic Nuclei, despite decades of research, are still not 
sufficiently understood. The basic theory of accretion disks gives us insight into physics 
involved in the process of accretion but is unable to explain the observed characteristic 
and time evolution of real accretion disks. Our understanding of the source of viscosity, 
spectral states, quasi-periodical oscillations, jets etc. is still far from satisfactory. New 
ideas and models are necessary to better understand one of the most powerful objects 
in the Universe. 

1.1 Origins of black hole physics 

As early as in the end of 18th century John Michell^ and Laplace noted that according 
to the Newton's theory of gravity and the corpuscular theory of light an object so 
massive and compact that even single photon is not able to escape to infinity may 
exist. However, this reasoning found only few supporters. 

In the early 1910's the general theory of relativity was established when Albert Ein- 
stein published a series of papers on gravitational field (Einstein, 1915). These works 
were followed in a short time by Karl Schwarzschild's solution of Einstein's equations in 
the vicinity of a spherically symmetric mass (SCHWARZSCHILD, 1916). Both scientists 
appreciated this discovery but at that time they were not aware of the fact that it is 
the complete solution for spacetime surrounding a non-rotating, uncharged black hole 
(BH). 

In 1930's Chandrasekhar discovered an upper mass limit for spherical stars built 
from degenerate matter (CHANDRASEKHAR, 1931). A few years later Baade & 
ZwiCKY (1934) proposed that supernovae represent the transition of normal stars 
to neutron stars - compact objects supported by the pressure of degenerate neutrons. 
This discovery was followed by Eddington's remark that existence of the Chandrasekhar 

^"If the semi-diameter of a sphere of the same density as the Sun were to exceed that of the Sun 
in the proportion of 500 to 1, a body faUing from an infinite height towards it would have acquired at 
its surface greater velocity than that of light, and consequently supposing light to be attracted by the 
same force in proportion to its vis inertiae (inertial mass), with other bodies, all light emitted from 
such a body would be made to return towards it by its own proper gravity." (Michell, 1784) 
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limit inevitably leads to a collapse of a sufficiently massive object to a BH (EddinG- 
TON, 1935). However, he considered such objects unlikely and tried to modify the 
degenerate equation of state to obtain stable, non-collapsed, solutions for any mass of 
a star. Such an approach was common among astronomers at that time. 

Oppenheimer & Snyder (1939) considered the collapse of a uniform, pressureless 
sphere of gas in the framework of the general theory of relativity. They analytically 
proved that an object resulting from the collapse is, at its particular stage, cut off from 
information exchange with the surroundings. For the ffist time an astrophysical process 
leading to the creation of a BH was proposed. 

Due to the lack of further research BHs remained outside the interest of most physi- 
cists until 1960's when Kerr (1963) solved Einstein's equations in the vicinity of a 
rotating, uncharged compact object. This solution was generalized by NEWMAN ET 
AL. (1965) to charged objects. Today we know that the Kerr-Newman geometry per- 
fectly describes the gravitational and electromagnetic fields of a stationary, solitary 
BHs. However, the importance of these discoveries were appreciated only a dozen or so 
years later. 

Since that time a large number of papers devoted to different aspects of the BH 
theory have been published. The research was motivated by new discoveries: quasars 
(Schmidt, 1963), pulsars (Hewish et al., 1968) and X-ray binaries (Oda et al., 
1971). The observations of the Cygnus X-1 microquasar became ffist evidence for the 
existence of BHs in the Universe. 

1.2 Black holes in the Universe 

Astrophysical BHs originate mainly through the process of a gravitational collapse. It 
takes place when the gravity force acting on a massive object is not balanced by other 
forces e.g., gas or radiation pressure. Gravitational collapse has been common in the 
Universe. Such a process led to the creation of separated clouds of matter after the Big 
Bang, which later transformed into clusters of galaxies and ultimately into galaxies, 
stellar clusters, clusters and planets. 

Gravitational collapse is especially important in the stellar evolution. Protostellar 
clouds of matter collapse leading to the increase of the pressure and temperature, and, 
as a result, to the ignition of the nucleosynthesis in their centers. Stars collapse also in 
the end of their evolution once their internal sources of energy have run out. If the core 
contraction is stopped by the pressure of degenerate electrons, a white dwarf originates. 
This scenario is possible only for low mass stars. If the Zero Age Main Sequence (ZAMS) 
mass of a star exceeds ~ 8Mq then the core exceeds the maximal mass of a white dwarf 
[IAMq, the Chandrasekhar mass) and continues its collapse. Only the pressure of 
degenerate neutrons can stop this process leading to a neutron star remnant. However, 
the cores of most massive stars (> 20Mq) exceed the Tolman-Oppenheimer-Volkoff 
mass (Oppenheimer & Volkov, 1939) and even degenerated neutrons are not able 
to prevent further collapse. Such an objects continues its contraction and, once it has 
fallen inside its own event horizon, becomes a BH. 

The scenario described in the previous paragraph leads to the creation of stellar- 
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mass BHs. The resulting mass of the remnant depends mostly on the initial (ZAMS) 
mass of a star and its metallicity. According to the current evolutionary models (e.g., 
Hurley ET AL. (2002)), masses of BHs resulting from evolution of single stars in 
the Galaxy range between 4 and 12Mq depending on the initial stellar mass. In low- 
metallicity environments, e.g., globular clusters, their masses may be much higher, 
reaching 25Mq, due to decreased efficiency of stellar winds. 

During the supernova explosion, BHs resulting from the stellar evolution may re- 
main in binary systems if the remnant does not obtain a kick sufficiently strong to 
disrupt a binary system or if the BH originates as a result of a direct collapse of the 
star. If the binary is close enough, the BH companion may fill its Roche lobe and 
start transferring mass onto the compact object. It may also happen that the compan- 
ion's envelope embraces the BH. During such an event, called the common envelope, 
the system circularizes or, if there is not enough orbital energy to eject the envelope, 
merges giving rise to a more massive single BH or a Thorne-Zytkow object (Thorne 
& Zytkow, 1977). 

As it will be pointed out in detail in this thesis, matter falling onto a compact 
object loses its angular momentum in an accretion disk and emits radiation. In case 
of stellar mass BHs this emission falls into X-rays and is amenable to observations by 
X-ray satellites (e.g., XMM-Newton or RXTE). Detailed observations of light curves, 
together with spectroscopic research, provided estimates of BH masses in a number of 
systems (Table 1.1). In case of these objects, the lower mass estimates are higher than 
2.5Mq, which makes them likely BH candidates. So far, these objects have been the 
only galactic BHs known. 

BHs of masses exceeding IO^Mq are called supermassive and form another class of 
astrophysical BHs. They are observed indirectly in nuclei of many galaxies, in particular 
in the center of the Milky Way. It is assumed that the evolution of most, if not all, 
galaxies leads to the creation of a very massive BH in their nuclei. 

The existence of supermassive BHs in active galactic nuclei (AGN) is supported by 
observations of visual/UV spectra of accretion disks which may be explained only under 
the assumption that they surround supermassive compact objects. In case of other than 
AGN types of galaxies we base on the Doppler measurements of water masers movement 
in the vicinity of nuclei. Their velocities reflect the existence of massive, dark objects in 
nuclei of galaxies. According to our knowledge, only BHs are compatible with provided 
characteristics. Finally, the observed orbital motion of stars in the vicinity of the 
Milky Way center (Sgr A*) suggests the supermassive BH of mass 4.1 ih 0.6 ■ IO^Mq 

(DOELEMAN ET AL., 2008). 

BHs of masses of the order of hundreds or thousands of solar masses form one more, 
so far only hypothetical, class of astrophysical compact objects - intermediate mass BHs. 
Observed Ultra-luminous X-ray Sources (ULXs) in nearby galaxies seem to suggest that 
BHs of such masses exist. However, the observed high luminosities may be explained 
by different processes e.g., by super-critical accretion onto stellar mass BHs. Doppler 
measurements of the velocities of objects gravitationally bound to intermediate mass 
BHs would provide evidence for their existence. However, so far no such observation has 
been obtained. It is also not clear how intermediate mass BHs would originate. They 
are too massive to result from the evolution of a single star of Population I or II. The 
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regions where they are found differ significantly from galactic nuclei so the processes 
giving rise to supermassive BHs cannot be involved. Intermediate mass BHs could be 
a product of hierarchical mergers of massive, unevolved stars into very massive objects 
collapsing into a BH. 

There is one more class of BHs that still remains hypothetical — mini BHs, with 
masses smaller than any BH resulting from a stellar collapse. It is possible that such 
quantum primordial BHs were created in the high-density environment of the early 
Universe (or big bang), or possibly through subsequent phase transitions. They are ex- 
pected to emit weak Hawking radiation, but direct detection of them will be impossible 
in the nearest future. 



Table 1.1: Galactic black hole candidates 



Name 


Other names 


Spectral type 
(0) 


Porh 


Mbh/M© 


Msec 


Refs. 


Cyg X-1 


HDE 226868 (0), 


09.7 lab 


— -J ^ 

5'^6 


13.5^29 


40 ± IOM0 


(1) 




V1537 Cyg (0) 












LMC X-1 




07-9 III 


3'^91 


10.91 ± 1.41 


31.79 ±3.48 


(2,3) 


LMC X-3 




B3 V 


l'^70 


6^9 


[0.5 ^ 0.91] 


(4) 


SS 433 


V1343 Aql (0) 


A3-7 I 


13'^1 


2.9 ±0.7 


1O.9±3.1M0 


(5) 


GX 339-4 


V821 Ara (0) 


F8-G2 III 


1*^756 


> 6 


[< 0.08] 


(6) 


GRO J0422+32 


V518 Per (0), 


M2 V 


5''09 


3.97 ±0.95 


[0.1161°;°??] 


(7) 




XRN Per 1992 












A 0620-00 


V616 Mon (0), 


K4 V 


7^75 


11 ±2 


[0.074 ± 0.006] 


(8) 




Mon X-1, 














XN Mon 1975 












2S 0921-630 


V395 Cor (0) 


KO III 


gd 


i.JU„o.25 • '^■-'-QA 


[1.32 ±0.37] 


(9,10) 


GRS 1009-45 


MM Vel (0), 


K8 V 




~ 8.5 


~ 0.5 


(11) 




XN Vel 1993 












XTE J1118+480 


KV Uma (0) 


K7-M0 V 




8.53 ±0.60 


0.09 ^ 0.5Mq 


(12) 


GS 1124-683 


GU Mus (0), 


KO-5 V 


lOM 


7.0 ±0.6 


[0.128 ±0.04] 


(13) 




XN Mus 1991 












GS 1354-645 


BW Cir (0) 


GO-5 III 


2'^54 


> 7.83 ±0.50 


1.07^ 2.11Mq 


(14) 


4U 1543-47 


IL Lup (0) 


A2 V 


l'^12 


8.4^10.4 


[0.25^0.31] 


(15) 


XTE J1550-564 


V381 Nor (0) 


G8 IV-K4 III 


l'^55 


9.7^11.6 


[0.09^0.15] 


(16) 


XTE J1650-500 




K4V 


0'^3205 


2.73^7.3 




(17) 


XTE J1652-453 








< 30 




(18) 


GRO J1655-40 


V1033 Sco (0), 


F3-6 IV 


2'^62 


6.3 ±0.3 


[0.40 ± 0.03] 


(19) 




XN Sco 1994 












MAXI J1659-152 






2Ml 


20 ±3 




(20) 



Table 1.1 - continued from previous page 



IN dllic 




(0) 


P , 

orb 


/ A/f 




iveio. 


n 1 ( uo-zou 


vziU( v^pn l^vJJ 


ivO V 


1 9'*p; 


7 • 7 Q 


[U.U14_oo^2J 


[Ai) 


XTE J1817-330 








< 6.0tt° 




(22) 


"VTTT' T1 81 Q OKA 


V ^tO^ti Ogl ^^^J J 


RQ TTT 
oy ill 


z oi ( 


R 8 • 7 /I 

O.O — ( .41: 


rn /1 9 • n /I 

[U.^tZ — U.^tOJ 


[Z6) 


XTE J1856+053 








< 4.2 




(24) 


XTE J1859+226 


V404 Vul (0) 


~G5 


9'^16 


7.6 + 12 




(25) 


GRS 1915+105 


V1487 Aql (0), 
XN Aql 1992 


K-M III 


33'^5 


10 + 18 


> 1.2M0 


(26) 


GRS 2000+251 


QZ Vul (0), 
XN Vul 1988 


K5 V 


8'^3 


5.5 + 8.8 


0.16 + 0.47 


(27) 


GRS 2023+338 


V404 Cyg (0), 
XN Cyg 1989 


KO IV 


Qd4Q 


10 + 13.4 


[0.060 ± 0.005] 


(28) 



Only objects with existing estimates of the BH mass are listed. 

(O) - optical counterpart, Porb - orbital period, Mbh - black hole mass. 

Msec - secondary mass or mass ratio (Msqc/Mbh) (in square brackets). 

References: (1) ZiOLKOWSKi (2004), (2) CowLEY ET AL. (1995), (3) Orosz et AL. (2009), (4) SORIA ET AL. (2001), 
(5) HiLLWiG ET AL. (2004), (6) Hynes et AL. (2003), (7) Gelino & Harrison (2003), (8) Shahbaz et al. (1994), 
(9) Shahbaz et al. (2004), (10) Jonker et al. (2005), (11) Macias et al. (2011), (12) Gelino et al. (2008), 
(13) Takizawa et al. (1997), (14) Casares et al. (2004), (15) Orosz et al. (1998), (16) Orosz et al. (2002), 
(17) Orosz et al. (2004), (18) Hiemstra et al. (2011), (19) Orosz & Bailyn (1997), (20) Kuulkers et al. (2011), 
(21) Remillard et al. (1996), (22) Sala et al. (2007), (23) Hjellming et al. (2000), (24) Sala et al. (2008), 
(25) Orosz (2003), (26) Greiner et al. (2001), (27) Ioannou et al. (2004), (28) Shahbaz et al. (1994). 
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1.3 Accretion onto compact objects 

In astrophysics the term accretion is used to describe the process of acquiring matter, 
typically gaseous, by a massive body. As a result of the infall, the gravitational energy 
is extracted. Its amount is roughly proportional to the ratio of the mass and radius 
of the attracting object: M/R. It is clear that the more massive and compact the 
object is, the larger amount of energy may be released. Accretion onto compact objects 
(BHs, neutron stars or white dwarfs) is therefore an extremely important process in the 
astrophysical context, as it transforms the gravitational energy into radiation in many 
of the most luminous classes of objects in the Universe. 

The spherically symmetric (zero angular momentum) accretion on a gravitating 
body is the simplest case one can imagine. It was first studied by BONDI (1952) and 
is called the Bondi accretion. It was shown that the solution describing such a process 
is unique and that the flow becomes transonic crossing the critical (sonic) radius. It is 
currently accepted that Bondi-type of accretion hardly occurs in the Universe because 
all potential sources of accreting matter (e.g., gas clouds or stars) have non-zero angular 
momentum. 

As an example let us consider a molecular cloud collapsing under the gravity of 
the central compact object. As it becomes denser, the random gas motions that were 
originally present in the cloud, average out in favor of the direction of the nebula's net 
angular momentum. Conservation of angular momentum causes the angular velocities 
to increase as the nebula becomes smaller. This rotation causes the cloud to flatten out 
and take the form of a disk. 

Formation of a disk is required to support the accretion rate. A free particle with 
a given angular momentum will rotate around the central object in an eccentric orbit. 
As long as angular momentum is conserved, the test particle will not fall onto the 
gravitating body nor release its potential energy. In such a case there is no accretion 
and no energy is released. To support the accretion and form a luminous accretion disk 
the angular momentum must be transported outward enabling particles to flow inward. 
Such a process occurs due to the action of viscous torques within the disk. It may 
be shown that laminar viscosity cannot provide sufficiently high magnitude of torque 
which would lead to the observed disk luminosities. Turbulent motions are required 
to explain the angular momentum transport in accretion disks. The exact mechanism 
leading to turbulence is discussed in Section 2.2. 

Accretion disks appear in many astrophysical objects. Protostellar clouds of molecu- 
lar gas create protoplanetary disks surrounding newly formed stars. Stars in close com- 
pact binary systems transfer mass through accretion disks when they fill their Roche 
lobe or due to strong stellar winds. Supermassive BHs in AGNs are believed to accrete 
cold matter from their vicinity. According to the collapsar model, short 7-ray bursts 
accrete the ejecta at super-critical rates shortly after the outbursts. Other examples of 
accretion disks in the Universe may also be given. Accretion appears to be involved in 
most of the high-luminosity phenomena in the Universe. 
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1.4 Spectral states 

The timescale for propagation of matter in accretion disks (the drift timescale) is much 
longer for disks around supermassive BHs than for disks in microquasars. For the 
latter, the changing mass supply from the companion may cause variability on the 
timescale of days. All of the observed microquasars exhibit this kind of variability 
proving that the process of accretion is not, in general, stationary. In the top panel 
of Fig. 1.1 we present, after DUNN ET AL. (2008), light curves in the hard and soft 
bands of microquasar GX 339-4 during the period March 2002 — July 2003. Both 
fluxes vary over a few orders of magnitude. They are, however, anticorrelated which 
means that the shape of the spectrum has to change significantly. The bottom panel 
presents luminosity vs hardness ratio. The microquasar exhibits the behavior typical 
for all objects of this kind: a given cycle starts with a hard spectrum of low luminosity 
(bottom right corner) which is followed by rapid increase of luminosity and softening 
of the spectrum leading to quasi-stationary state characterized by soft spectrum and 
high luminosity (top left). Close to the end of the cycle, the microquasar decreases 
its luminosity preserving, however, soft spectrum (left branch). Finally, the spectrum 
hardens again, but this time the change of hardness takes place at lower luminosity 
than during the beginning of the cycle. 

Two dominant spectral states can be distinguished: high/soft (soft spectrum and 
high luminosity) and low/hard (hard spectrum, low luminosity). Fig. 1.2 presents 
schematic pictures of the spectral energy distributions corresponding to these states on 
left and right panels, respectively. During the high/soft state the spectrum is domi- 
nated by the soft component which may be well modeled by superposition of black body 
radiation with different effective temperatures and magnitudes (multicolor black body) 
which corresponds to the emission predicted for an accretion disk (see e.g., Eq. 3.18). 
The high energy tail comes most probably from photons up-scattered by hot electrons 
above the disk. The low/hard state is, on the contrary, dominated by the hard com- 
ponent, while the disk radiation contribute only to small fraction of the total emitted 
energy. 

The disk structure behind the high/soft state is relatively easy to resolve — the 
spectrum can be explained by an optically thick accretion disk extending down to 
the BH horizon, which emits the viscosity-generated heat as multicolor black body 
radiation. The low/hard state, on the contrary, requires much more complicated disk 
structure. To obtain the observed relation between the hard and soft components one 
has to assume that the standard, optically thick accretion disk terminates far away 
from the BH, at radii of the order of hundreds of the gravitational radius. Inside 
this transition radius the flow has to be radiatively inefficient and hot so that it can 
up-scatter photons originating in the outer disk regions. 

A self-consistent model explaining the characteristic features of the low/hard spec- 
tral state has not been developed yet. The transition from the optically thick to thin 
disk may be a result of disk evaporation due to conductive heating from the corona 
(e.g., Meyer ET AL. (2000)). The high/soft state, on the contrary, is well understood 
and consistent models of the underlying optically thick disks have been in use for many 
years. In this dissertation we address this state of accretion disks only. 
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Outburst -J 




Figure 1.1: Soft (3 - lOkeV) and hard (6 - lOkeV) flux components for GX 339-4 
during one cycle of observations (top panel). Luminosity plotted against hardness ratio 
for the same set of observations (bottom panel). Colors denote particular observations. 
Figures after Dunn et al. (2008). 
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energy energy 



Figure 1.2: Schematic plots presenting general features of microquasar spectra during 
high/soft (left) and low/hard (right) spectral states. The red curve corresponds to the 
disk component, the purple to the comptonized component, while the green curve fits 
the iron line at 6.4 keV. 



Chapter 2 

Accretion disk theory 



2.1 Conservation laws 

The flow of matter accreting onto a compact object tlirougli an accretion disk may be 
described by tlie equations of hydro- (or, more general, magnet ohydro-) dynamics as 
the typical mean free path of collisions of microscopic particles in such medium is much 
shorter than the characteristic macroscopic length of the disk (e.g., its thickness). 
The first fundamental equation describing fluid dynamics is, 

V{ri = 0, (2.1) 

where the hydrodynamical stress-energy tensor has the general form^, 

Tl = u'ukip + e) + p5l - tl + u'qk + Ukq\ (2.4) 

where is the matter four-velocity, p is the total pressure, e is gas internal energy, 61 
is the Kronecker delta, is the radiative energy flux and the viscous stress tensor t], is 
given by, 

4 = i^P^l (2.5) 

(z/ is the kinematic viscosity coefficient and al is the shear tensor of the velocity field). 

The second fundamental equation is that of the conservation of the baryonic number 
of particles, 

Viipu') = 0, (2.6) 

where p is the gas density. 

^In the case of magnetohydrodynamics, when the viscosity resuhs directly from magnetorotational 
instabihty, as discussed in Section 2.2.2, the stress-energy tensor is, 

n = n'ukip + e + b^) + {p + }? 12)^1 - b%k + u'qu + Ukq\ (2.2) 

The contravariant fluid- frame magnetic field is given by 6% and is related to the lab- frame three-field 
by, 

V ^ B^Kl/u' (2.3) 

where h\. = u^Uk -|- (5^ is a projection tensor. 

Throughout this dissertation we use the ( — h -t--f) signature. 
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These general equations of hydrodynamics are used to construct all relativistic, 
hydrodynamical models of accretion disks. In Appendix B we give a detailed deriva- 
tion of the non- stationary equations describing advective, optically thick, vertically- 
integrated, a-disks. In this section we introduce and discuss, basing mostly on Kato 
ET AL. (2008), the equations governing the stationary, axisymmetric accretion flow in 
the non-relativistic limit. 

Under this approximation the fluid dynamics is described by the Navier-Stokes equa- 
tions: 

|^ + V-(pn) = 0, (2.7) 
p(^ + (m-V)m) = -pV<l> - Vj9 + Vt, (2.8) 



where u is the fluid velocity, $ is the gravitational potential and t is the viscous stress 
tensor. 

The surface layers of accretion disks usually do not extend far from the equatorial 
plane. It is therefore possible to use the cylindrical system of coordinates (r, (j), z) . The 
equation of continuity takes the form. 

We neglect both time and azimuthal derivatives and, following the common approach, 
we integrate this equation vertically between the disk vertical boundaries at z = ±h, 

1 d d 



which is equivalent to, 

ld_ 

r dr 



/ purdz + puz\th = ^- (2-11) 



As the density p is very low at the disk boundary, the second is negligible. Therefore, 
we may write, 

|-(rSV^) = 0, (2.12) 



where 

fh 



pdz, (2.13) 

h 



1 



s 



V = - purdz, (2.14) 



h 



are the surface density and density-weighted radial velocity, respectively. Finally, inte- 
grating Eq. 2.12 we may introduce the integration constant M, 

M = -27rrSV, (2.15) 

which is the accretion rate of gas. This equation describes the conservation of mass. 



2.2 Source of viscosity 
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Taking the (j) component of Eq. 2.8 we get, under the same assumptions: 

pur-^nr^ = i|-(r\^), (2.16) 
or r or 

where Q is angular velocity of the fluid and is the only non-negligible component of 
the viscous stress tensor. Introducing the angular momentum per unit mass C = Qr^ 
and integrating over z and r we obtain the angular momentum conservation law: 

^iC-CJ = -2hr\^. (2.17) 
Zn 

We have assumed that there is no torque at the inner boundary (i.e., BH horizon) so 
that the integration constant Ci^ can be interpreted as the angular momentum of the 
matter at the inner edge. Eq. 2.17 reflects the conservation of angular momentum 
which is transported outward by viscous stresses. 

Energy is another conserved quantity. Three mechanisms of heating/cooling are 
involved for advective accretion disks: viscosity, advection and radiation. The energy 
balance may be formulated in the following way: 

gadv _ gvis _ grad (2.18) 

which reflects the fact that the heating or cooling rate by advection {Q^'^^) must be 
equal to the difference between the amount of heat generated by viscosity (Q"*^) and 
the radiative cooling rate {Q'^^'^). 



2.2 Source of viscosity 

The angular momentum in accretion disks is transported by viscosity. In this section 
we will discuss this process in general, describe the ap formalism and discuss the mag- 
netorotational instability which is considered the most likely source of turbulence in 
disks. 

Let us consider parallel shear flow (left panel of Fig. 2.1). As the upper layer moves 
with greater velocity the momentum transfer occurs from the top to bottom. As a 
result the upper layer will be decelerated while the bottom accelerated. Viscosity tends 
to reduce shear and produce uniform layers. 

The viscous force exerted in the horizontal {x) direction on the unit surface of the 
interface plane is, 

t.y = V^, (2.19) 

where y is the normal to the interface, t] = pu is the coefficient of dynamical viscosity 
and dvx/dy is the velocity gradient. 

In the case of accretion disks we do not deal with parallel layers. Instead, we have 
differentially rotating annuli (right panel of Fig. 2.1). The rotational velocity increases 
in Keplerian flows inward. Therefore, the viscosity will transfer angular momentum 
from the inner annulus to the outer. In such a way, matter loses angular momentum 
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Figure 2.1: Parallel shear flow (left) and differentially rotating flow (right panel). Figure 
after Kato ET AL. (2008). 



and moves inward. If the disk was cut off at a certain outer radius and no external 
force acted on the disk, the outermost annuli would receive all the transported angular 
momentum and move outward to infinity. 

The viscous force per unit area exerted in the azimuthal direction of the interface 
plane between differentially rotating annuli is, 

i,, = „(^-l>)=„.f?. (2.20) 

\ or r J or 

2.2.1 ap prescription 

In 1970s the understanding of viscous processes which take place in accretion disks 
was very poor. Nevertheless, Shakura & SUNYAEV (1973) proposed an extremely 
elegant and plausible formalism which has been used effectively for many years, making 
progress in accretion disk theory possible. They assumed correctly that the source of 
viscosity in accretion disk is connected in some way to turbulence in gas-dynamical 
flow. The nature of this turbulent motion was not known at that time. They assumed 
a simple expression for kinematic viscosity coefficient: 

U Wturb/turb, (2.21) 

where Vturh and /turb stand for typical turbulent motion velocity and length scale, re- 
spectively. Assuming that the velocity of turbulent elements cannot exceed the speed 
of sound Cs and their typical size cannot be larger that the disk thickness h one gets 

1/ < Csh. (2.22) 

Using the vertical equilibrium condition (Eq. 3.13), this expression may be put in the 
following form 



The stress is given by Eq. 2.20: 



tr4> = W-Q^ ~ -pvVt. (2.24) 



2.2 Source of viscosity 
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Taking Eq. 2.23 into account one gets 

\tr^\ < p^n = p. (2.25) 

Therefore, a non-dimensional viscosity parameter a satisfying the condition a < 1, may 
be introduced: 

tr<j> = —ap- (2.26) 

This expression for viscosity, proposed for the first time by Shakura & SUNYAEV 
(1973), is called the ap or the a prescription and has been widely used for many 
years in accretion disk theory. The most modern magnetohydrodynamical (MHD) 
simulations (e.g., HiROSE ET AL. (2009b)) support application of the a formalism to 
mean accretion flow dynamics. 



2.2.2 Magnetorotational instability 

For almost twenty years the source of the turbulent motion responsible for viscosity in 
accretion disks remained unknown. Only in 1991 BalbuS & Hawley (1991) identified 
the process which might lead to turbulent motion in accretion disks. They noted that 
in astrophysical disks a magnetorotational instability (MRI) may set in, amplifying the 
initial weak magnetic field and leading to the required angular momentum transport. 
Since that time a number of numerical methods have been developed to verify this 
possibility. They gave positive results and therefore magnetorotational instability is 
currently thought to be the most likely explanation of turbulent motion in accretion 
disks. However, it was suggested by Umurhan et AL. (2007) on the basis of perturba- 
tive weakly nonlinear analysis that the transport due to the nonlinearly developed MRI 
may be very small in experimental setups with the Prandtl number Pm <^ 1- Several 
authors studied an alternative picture of turbulence arousal in non-magnetized astro- 
physical disks, related to their non-axisymmetric hydrodynamical disturbances, which 
can, in principle, precipitate secondary instabilities (Rebusco ET AL., 2009). 

MRI was first noticed by Velikhov (1959) who studied stability of an ideally 
conducting liquid flowing between cylinders rotating in a magnetic field. His results 
were generalized by Chandrasekhar (1960). Only in 1990's astrophysicists realized 
its paramount importance for triggering accretion in astrophysical objects. 

The general idea of MRI may be illustrated in the following way. Let us consider 
a disk of fluid rotating with Keplerian angular velocities in a weak axially symmetric 
magnetic field. The magnetic tension acting on two fluid elements may be imagined 
as a spring connecting them. As the inner element orbits with higher angular velocity 
than the outer one, the spring stretches. In such a way the magnetic tension slows down 
the inner element and forces the other one to speed up. As a result, the inner loses its 
angular momentum while the angular momentum of the outer element increases moving 
the elements apart. In this way the angular momentum is transported outward making 
accretion possible. The crucial condition for arising of MRI is that the angular velocity 
of a magnetized fluid decreases with radius: 



(2.27) 
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Figure 2.2: Time dependence of the ratio ot of the box-averaged r</) stress to various 
box-averaged pressures in one of the simulations from HiROSE ET AL. (2009a). From 
top to bottom, the stress is scaled with gas pressure, the geometric mean of gas and 
radiation pressures, and with total (gas plus radiation) pressure, respectively. Figure 
after HiROSE ET AL. (2009a). 



which is satisfied in the astrophysical context, e.g., for the Keplerian profile of rotation. 

Most of research devoted to MRI is based on numerical simulations solving the 
equations of magnetohydrodynamics. They may be solved either globally in three- 
dimensions (e.g., Penna ET AL. (2010)) or locally in the shearing-box approximation 
(e.g., HiROSE ET AL. (2009a)). So far, the radiative transfer has been neglected as 
it immensely increases the required computational power. Instead, different, purely 
numerical, methods of cooling are assumed. Only recently, first MHD simulations 
accounting for radiation, however axially symmetric, were performed (Ohsuga et 
AL., 2009). 

The crucial question from the point of view of hydrodynamical models of accre- 
tion disks is whether the assumption made by Shakura & Sunyaev remains valid facing 
the MHD physics. Fig. 2.2 shows time-dependence of the value of a (defined as the 
ratio of the box-integrated stress to the different box-integrated pressure quantities) for 
various stress prescriptions for simulation 1112a from HiROSE ET AL. (2009a). The 
short-time variability reflects the fluctuations in the r</) stress component. The total 
(gas plus radiation) case has the least variation of the a profile (by a factor of 4), while 
the largest (by a factor of 10) is for the gas pressure only case. Therefore, the stress 
component never follows the pressure precisely, but one may say that the prescrip- 
tion (with "p being the total pressure) quite well describes the mean dynamics of the 
flow. 



2.3 Efficiency of accretion 
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2.3 Efficiency of accretion 

Let us consider the energetic balance of the accreted fluid between two radii: and 
Tout- The amount of generated energy (luminosity) between these two radii can be 
estimated from 

^nn^ro.. = Mer^ut + Crfi)w - Mer,^ - (rfi).,„, 

where T denotes the torque at a given radius, Tf2 is the outward energy transport due 
to viscosity and e is the specific energy. If the disk extends to large radii, we may put 
fir t ~ and Cr » ~ 1. Hence, 

' out ' out } 

Lr,^ = M (1 - e,.„ - xi^rj = VM, (2.28) 

where x = T/MC is the ratio of the viscous torque to the advective flux of angular 
momentum (see Figs. 4.18 and 4.19). 

If the torque at the inner edge of the disk is negligible {x ^ 1)) the efficiency of 
accretion fj depends mainly on the specific energy at the inner edge, e^j^, given by^ 



2M 

1 - — . 2.29 

3ri, 



m 



Assuming x = we get, 

fj = l- e,,„. (2.30) 

The farther away the inner edge from the innermost stable circular orbit^ (ISCO, Ap- 
pendix A. 3) is (and the closer to the BH), the smaller the efficiency. 
For accretion disks terminating at ISCO, we have. 



l-,/f = 0.057, a, = 0, 

(2.31) 

1 - V I = 0-420, a, = 1. 

Fig. 2.3 presents the dependence of the efficiency of accretion on the BH spin parameter 
a*. For non-rotating BHs it equals rj = 0.057. The efficiency rapidly increases for BH 
spin values close to unity — even small difference in angular momentum (e.g., a* = 
0.999 vs 0.9999) implies significant difference in the efficiency of extracting gravitational 
energy from the flow. 



^Throughout this dissertation we use the geometrical {G = c = 1) units. In these units mass is 
measured in meters, M = = GM/c^ and Mq k, 1477 m. 

^The innermost stable circular orbit (ISCO) is often called the marginally stable orbit. 
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Figure 2.3: Efficiency of accretion {fj = L/Mc^) through a standard thin disk termi- 
nating at ISCO for a given value of BH spin. 



Chapter 3 

Classical models of accretion disks 



Figure 3.1 shows the locations of analytic and semi-analytic models of hydrodynamical 
accretion disks^, in the parameter space described by the vertical optical depth r, 
dimensionless vertical thickness and dimensionless accretion rate m = M/MEdd (defined 
in Eq. 4.35). In the following sections we describe in detail Polish doughnuts and 
radiatively efficient disks (Shakura & Sunyaev, 1973; NoviKOV & Thorne, 1973). 
Slim disks are introduced and discussed in Chapters 4 — 8. 

3.1 Polish doughnuts 

One of the simplest and most useful models of geometrically thick accretion disks may 
be obtained adopting the following assumptions (Abramowicz ET AL., 1978): the 
matter is distributed uniformly and axisymmetrically, it orbits the BH along circular 
trajectories with given angular momentum L = u^{r, 6) and the stress energy tensor 
is given by: 

Ti = u'u.ip + e) + 5lp (3.1) 

where u^^ is the matter four-velocity, p the total pressure, e the total energy density and 
5^ is the Kronecker delta. 

The four velocity of a particle on a circular orbit has the following form, 

u' = A [t]' + nC) , (3.2) 

where r]"^ and ^* are the Killing vectors associated with the time and axis symmetry, 
respectively (Appendix A.l) and Q = u'^/u^ is the angular velocity. The normalization 
constant A comes from the u^Ui = — 1 condition, 

= g^^ + 2ngt^ + n^g^^. (3.3) 

From the general equation VjT^ = it follows, 

= V.T^ = Vi {u'ukip + e) + 61p) = UfcV. {{p + ey) + u'(p + e)ViUk + VkP- (3.4) 

For detailed reviews of these solutions SGG, http: //www, scholarpedia . org/ article/ Accretion_discs OT Kato et al. 

(2008). 
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optical depth at 20 M 



Figure 3.1: This figure illustrates a few of the most well-known analytic and semi- 
analytic solutions of the stationary black hole accretion disks. Their location in the 
parameter space approximately corresponds to viscosity a = 0.1 and radius r = 20 M. 

The first term vanishes as it reflects the conservation of number of particles. Thus, we 
have, 

__Vfcp_ ^ u'ViUk = ttk- (3.5) 
p + e 

The acceleration term may be expressed as follows, 

ak = Au' (V,r/fc + nV.ik + ^fcV.fi) = Au' {-VkVi - ^^k^i + ^k^i^) = (3-6) 
= -Au'Vk{v^ + ^ii) + Au' i^iVk^ + ^kViQ) = 

= -^^'Vfc {gtt + 2Qgt^ + Q^g^^) + Au' {^iVk^ + ^Vil^) = 

= -^A^ {Vkgtt + '2nVkgt4> + ^^^^kg^) - \a^ {'^gt<iNk^ + 2g^^nVk^) + 

+ Au' i^iVk^ + ^k^iQ) = ~A^ [Vkgtt + ^QVkgt^ + n^Vkg^^) + A^^ku'ViQ. 
The last term vanishes as f2 = Q{r). In the above derivation we have used the Killing 
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equation (VfcT^j + Virjk = 0). Finally, we get 

Vfcp 1 



p + e 2 



{Vk9tt + 2nVk9tr^ + n^Vk9^^) ■ (3.7) 



By rewriting the r and 6 components of this equation and dividing them side by side we 
may get the following expression for the shape of the equipressure surfaces (VkP = 0), 

dr dgQu + 2Q dgg^ + deg^^ ' 

Knowing the angular velocity (or the specific angular momentum i = i{r, 9) = —u^/ut) 
one can integrate the function given above and obtain the exact profiles of equipressure 
surfaces. QiAN ET AL. (2009) have recently made a comprehensive study of their 
shapes. They parametrized the angular momentum distribution by introducing the 
following function, 

i(r,9) = I ^^K(risco) ^^^'"^ ^ ^ ^isco, (3.9) 

I ^{r = risco) sin2'^6' for r < risco, 

where risco is the radius of the marginally stable orbit. /3 — )■ corresponds to the limit 
of constant (in radius) angular momentum, while /3 = 1 and rj = 1 gives the Keplerian 
profile (^k('")) at the equatorial plane. The vertical distribution of angular momentum 
depends on the value of 7. Fig. 3.2 presents profiles of the equipressure surfaces in the 
Polish doughnut model for a non-rotating BH. Rows correspond to increasing value of 
P from (top) to 0.99 (bottom). Columns show 7 (0.0 ^ 0.9) dependence. The upper 
left panel shows a "standard" Polish doughnut. The lower right panel shows an almost 
Keplerian disk at the equatorial plane, surrounded by a very low angular momentum 
envelope. 

Fig. 3.3 presents the comparison of equipressure surfaces profiles with results of a 
modern MHD simulation (Fragile et AL., 2007a, b). It is clear that appropriate 
choice of parameters describing the analytical model may reproduce most of the rele- 
vant features of the numerical results, including the locations of the cusp and pressure 
maximum, as well as the vertical thickness of the disk. The discrepancy inside the cusp 
(the numerical solution maintains constant vertical height while the analytical equi- 
pressure surfaces diverge) is due to significant radial velocities in this region which are 
not consistent with the adopted assumption of circular motion of fluid. 



3.2 Radiatively efficient, optically thick disks 
3.2.1 Shakura & Sunyaev (1973) 

When discussing the history of models describing disk accretion onto a compact object 
one has to mention the outstanding, milestone work of Shakura & SUNYAEV (1973). 
The importance of this work was mainly due to innovative approach to viscosity which 
was described in details in Section 2.2.1. The authors introduced several assumptions 
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Figure 3.2: Equipressure surfaces in (r,z) coordinates for a non-rotating BH and 77 = 
1.085. Five rows correspond to /3 = (0.0), (0.1), (0.5), (0.9), (0.99) from the top to the 
bottom. Four columns correspond to 7 = (0.0), (0.1), (0.5), (0.9) from the left to the 
right. The black area in the bottom left corner of each plot denotes the BH. Figures 
after QiAN ET AL. (2009). 
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Figure 3.3: Comparison of pressure distributions between the analytic model {dark 
lines) and numerical simulations (colors). The results of MHD simulations (FRAGILE 
ET AL., 2007a, b) have been time-averaged over one orbital period at r = 50M. Upper 
panel: Schwarzschild black hole [a^ = 0); the analytic model parameters are rj = 1.085, 
P = 0.9, and 7 = 0.18. Lower panel: Kerr black hole (a* = 0.5); the analytic model 
parameters are r] = 1.079, (3 = 0.7, and 7 = 0.2. Figures after QiAN ET AL. (2009). 
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which are common for most of models describing disk accretion: (1) the companion 
star in a close binary has negligible gravitational influence on the disk; (2) the self- 
gravitation of the disk is not important; (3) disk is axisymmetric {d^ = 0) (4) disk 
is thin {h{r) ^ r) (5) disk is in stationary state {dt = 0), (6) particles move around 
the compact object on Keplerian circular orbits and (7) disk is radiatively efficient 
(radiative cooling is the only cooling mechanism). 

The physics of disk accretion under the assumptions given above is described by 
four formulae: conservation of rest-mass and angular momentum, vertical equilibrium 
of forces and energy balance. First two have already been derived in Section 2.1: 

(i) conservation of rest-mass (Eq. 2.15) 

M = -27rrSV, (3.10) 

where M stands for mass accretion rate (which is constant due to dt = assumption), E 
is surface density defined as S = J^^ p dz and V is the density- averaged radial velocity 
of matter. 

(ii) transport of angular momentum (Eq. 2.17) 

^{C - Ci,,) = -2hr\^, (3.11) 

where, assuming that there is no torque at the inner edge of the disk, the integration 
constant Ci^ can be interpreted as the angular momentum of the matter at the inner 
edge (£in = y/GMr^^ for the Newtonian potential), C = Qr"^ is the angular momentum 
and tj.^ is the only non-negligible component of the viscous stress tensor. 

(iii) vertical momentum conservation 

According to the assumption that the disk is in the stationary state, it reduces to a 
simple hydrostatic equilibrium condition. Writing down the balance of the vertical 
component of the gravitational force and the vertical pressure gradient we get: 

(3.12) 

pOZ j-^ J. 

Following the standard approach, we replace the differentials with finite differences 
{Ap ^ p{z = h) and Az h) obtaining the formula for disk semi-thickness h: 



where Cg is the sound speed in the midplane of the disk. 

For a Newtonian accretion disk, the t,.^ component of the viscous stress tensor is 
given by (compare Eq. 2.20), 

3 

Ucf, = --pu^, (3-14) 

where u is the kinematic viscosity coefficient which has not been defined yet. It can 
be shown (e.g., SHAPIRO & Teukolsky (1983)) that the heat is generated locally by 
viscosity at the rate: 

vis ^ M!^ (3^^5) 
up 
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Using Eq. 3.11 we can express trtf, in the following way: 

M 



(C-Cin). (3.16) 



Putting Eqs. 3.14 and 3.16 into Eq. 3.15 we obtain the total heat generation rate at 
radius r: 

0- = 2.,™ = |^l^fl-./li). (3.17) 

Assuming that the heat generated by viscosity is immediately radiated away we get the 
following expression for the outcoming flux of energy from one face of the disk (top or 
bottom): 

1 3M GM 

2 ovrr^ r 




= ( 1 - a/- ) ^ (3-18) 



which is the energy conservation law in the disk. It is of extreme importance to note 
that the expression for the outcoming flux does not depend on the viscosity prescription 
which we have not determined yet. This fact will have major consequences in estimat- 
ing BH spin basing on X-ray spectra (see the appropriate section or Shafee ET AL. 
(2006)). 

Let us now precise the treatment of radiation transport. The simplest thing one 
can do is to take the diffusive approximation (Section 6.1) and replace differentials with 
finite differences (as we did when formulating the vertical equilibrium): 

F{r,z) = -- — ^=^r^ (3-19) 

6Kp oz 3kL 

where k is the Rosseland-mean opacity coefficient and Tq denotes the temperature at 
the equatorial plane. In this thesis we often use the Kramer's formula: 

K = Kes + Kff = 0.34 + 6.4 X lO^V^""^-^ cmVg. (3.20) 

The total pressure of mixture of gas and radiation is given by, 

p = NkT + -aT^. (3.21) 

3 

Accretion disks are typically supported by the gas pressure with exception to the in- 
nermost regions where radiation pressure dominates. The formula given above assumes 
that the magnetic pressure is negligible. 

To close this set of equations, the (r, 0) component of the viscous stress tensor has 
to be specified. Shakura & Sunyaev (1973) introduced a non-dimensional viscosity 
parameter a satisfying the condition a < 1 (Section 2.2.1), so that, 

tr(p = —<yp- (3.22) 



This expression for viscosity is called the a prescription and has been widely used for 
many years in accretion disk theory. 
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The equations given above form a set which can be solved algebraically, e.g., for 
S(r), h{r) and T(r). It turns out that an accretion disk can be divided into three 
regions: (i) an outer region dominated by gas pressure and "free-free" absorptions, (ii) 
a middle region in which gas pressure and electron scatterings dominate and (iii) an 
inner region dominated by radiation pressure and electron scatterings (for low enough 
accretion rates, the middle and inner regions may not occur). Full formulae describing 
accretion disks in all three regimes may be found in the original work, as well as in, 
e.g., Kato et al. (2008). 

The seminal paper by Shakura & SUNYAEV (1973) laid the foundations of the 
theory of accretion disks. For the first time astrophysicists were able to investigate their 
structure. The ap prescription, introduced by the authors, has remained in wide use in 
astrophysics. 

In Fig. 3.4 we present a schematic picture of a radiatively efficient, geometrically 
thin accretion disk. 



3.2.2 Novikov & Thorne (1973) 

The obvious upgrade of the model given by Shakura & SUNYAEV (1973) was to write 
its equations in full general relativity. It was done a few months later by NOVIKOV 
& Thorne (1973). They considered accretion in the equatorial plane of a BH in 
the framework of the Kerr spacetime metric. The authors introduced a convenient 
formalism of splitting formulae into Newtonian limits times relativistic corrections. 
However, to retain compatibility with other relativistic models we describe in this paper, 
we will present NOVIKOV & Thorne (1973) disk solutions using formalism proposed 
in Abramowicz et al. (1996). 

Throughout this paragraph we will use the Kerr metric in cylindrical coordinates 
close to the equatorial plane (see Appendix A). We use geometrical units {c = G = 
1) and the following form of the stress-energy tensor of matter in the disk (compare 
Eq. 2.4): 

= pu'u'' + pg'^ - t"' + + u'q\ (3.23) 

where p is the rest mass density, f'' is the viscous stress tensor and is the radiative 
energy flux. All physical values in this paragraph, if not stated otherwise, are averaged 
in the same sense as in Section 2.1. We assume that the density of gas internal energy 
is much smaller than the density of gravitational binding energy. 

The rest mass conservation law requires that the covariant derivative of the product 
of rest mass density and four velocity disappears: 

V.(pn^) = ^(v^pn^),, = -irpu')^r = 0. (3.24) 

Integration over r introduces an integration constant which can be interpreted as the 
mass accretion rate: 

M = -27rSVAi/^ (3.25) 

where S = J^j^ p dz is disk surface density, A = — 2Mr + a^, and V = -^j^'u!' is now 
the fluid radial velocity as measured in the local rest frame. 




Figure 3.4: A scheme of the standard thin disk modeh Gas is accreted through a geometrically thin disk onto a BH. Viscosity 
results from turbulent motions. Orbits are nearly Keplerian. Disk terminates at the marginally stable orbit. Viscous heating is 
balanced by radiative cooling. Disk emits black body radiation with the effective temperature depending on radius. 
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The general form of the angular momentum conservation law is: 

V,{Tie) = 0, (3.26) 

where C,'' = S^^^ is the Kerr metric KiUing vector along (p direction (Appendix A.l). 
After some algebra we derive: 

-S^AV^^ - liL (r r\rA + 2FC = 0, (3.27) 



r dr r dr 



h 



where C = stands for the specific angular momentum per unit mass for circular 
orbit (please note that £ = —u^/uf is a constant of motion for perfect fluid) and F 
is the outcoming flux of energy from one face of the disk. It can be shown (Lasota, 
1994) that the component of the stress tensor is related to the comoving t^^ by the 
following relation: 

^1/2^1/2 

= ^^f^- (3.28) 

The Lorentz 7 factor is defined as: 

7-2 = 1- {Qrf, (3.29) 

where 

^ ^ ~ r3/2 + aMV2 (3.30) 

is the angular velocity of the Keplerian orbits in the Kerr metric. Let us denote the 
integrated shear stress J^l^ t^^ dz by T^^. Now, the angular momentum equation takes 
the following form: 

IeFAI/^^ --^( -A^'^^'I^^Tr^ + 2FC = 0. (3.31) 
r dr r dr yr J 

Similarly like in the Newtonian case, one can prove the following relation between shear 
and stress tensor components and the flux emerging from one side of the disk (NOVIKOV 
& Thorne, 1973): 

1 r^^ 1 



^.l^ j ^ t-^ dz = -af^Tr,, (3.32) 

where is the shear of the equatorial geodesic orbits which is equal to (Gammie & 
POPHAM, 1998): 

-fi = 2^^- (3-33) 
Putting the definition of flux from Eqs. 3.32 and 3.33 into Eq. 3.31 we finally obtain: 
d /Ai/2Ai/2^^ \ ^^2^^^^ Md£ 
dr \ r ^'^J r^ dr ^'^ 2tt dr 

The only unknown function is Tr<j){r). Therefore, the above equation can be directly 
integrated. Originally, NoviKOV & Thorne (1973) gave in their work a formula for 
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Trtf, involving integrals. However, the following paper of PAGE & Thorne (1974) gave 
explicit, algebraic expression for this quantity. Once we know the integrated shear 
stress, we can use Eq. 3.32 to calculate outcoming flux. To resolve radial dependencies 
of V, S, Tc etc. we must introduce formulae for viscosity, vertical structure and opacity. 

The authors followed Shakura & SUNYAEV (1973) and assumed that the comoving 
(r, 0) component of the viscous stress tensor can be expressed as (Section 2.2.1) 

tf^ = -Oip, (3.35) 

where a < 1 and p is the total pressure. This assumption leads to, 

Tr^ = -aP, (3.36) 

where P = J^l^ pdz is the vertically integrated pressure. 

The vertical hydrostatic equilibrium is constructed assuming that the pressure gradi- 
ent is balanced by the vertical component of the gravitational force which corresponds 
to the vertical epicyclic frequency Q±, 

Idp 2 M-H 

-— = -Q^z = -—z, 3.37 

poz r-^ C 

where "H and C are relativistic correction factor defined in Eq. 4.3. Once again, we 
replace differentials by finite differences to obtain: 

To close this set of equations we need to describe the energy balance. Similarly like in 
the previous section we assume that radiative transport is the only cooling mechanism. 
Introducing finite differences into the diffusive approximation we get: 

where r is the total optical depth. 

Taking all the equations given above together we get a set which can be solved in 
algebraic way, like in the Newtonian case. Similarly, we can distinguish three regions in 
the accretion disk depending on the dominating pressure and opacity sources. Explicit 
expressions for outgoing flux, surface density, disk semi-thickness, central temperature 
etc. are given in NOVIKOV & Thorne (1973) with additional formulae in PAGE & 
Thorne (1974). 
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Chapter 4 

Relativistic, stationary slim accretion 
disks 



Standard models of thin disks (Shakura & SUNYAEV, 1973; NOVIKOV & Thorne, 
1973) assume that all the heat generated at a given radius by viscosity is immediately 
radiated away, meaning that the thermal equilibrium is determined by the local balance 
of viscous heating and radiative cooling, = Q^'^'^. No other cooling mechanisms are 
considered in the standard model. This assumption is fairly consistent and correct 
only for sufficiency small accretion rates. For high accretion rates it fails as the simple 
argument below shows. 

The "advective" radial flux of heat, Q^*^^, is a consequence of the radial motion of 
hot matter. It can be estimated as (Kato et AL., 2008), 

g^'^-^-l\/p = -^^. (4.1) 

Assuming P = Pj-ad = 2H^Tq and taking the rate of radiative cooling, Q'^^'^ = 
we may estimate this expression in the following way, 

Q-v ^ JL.P , Q,.. /e^c = ^^Q-. (4.2) 

It is then obvious, that for M > {h increases with M), advective cooling pre- 

dicted by the Shakura & Sunyaev (1973) model is larger than the radiative cooling 
Q'^'^'^, and that the standard thin accretion model cannot be adequate at such accretion 
rates. Detailed calculations show that advection starts to modify the disk structure at 
accretion rates corresponding to L ^ O.GLEdd (Section 4.4). 

AbramowiCZ et AL. (1988) constructed optically thick "slim" accretion disk model 
which accounts for the advective cooling and does not impose the Keplerian angular 
momentum. Therefore, it may be applied to disks at any accretion rate. 

Advection appears also at much smaller (sub-Eddington) accretion rates, but for 
optically very thin disks. These disks, today known as advection-dominated accretion 
flows (ADAFs), were propheted in a paper by ICHIMARU (1977) (see also Rees et 
AL. (1982)). Only after their rediscovery in the mid 1990s by Narayan & Yl (1994, 
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1995) and by Abramowicz ET AL. (1995), ADAFs started to be intensely studied 
by many authors (see reviews in Lasota ( 1999a, b) and Narayan & McClintock 
(2008)). 

In this doctoral thesis, we study optically thick slim disks only, and this Chapter is 
devoted to the discussion of their relativistic, stationary model. 

4.1 Equations 
Assumptions 

We assume an axisymmetric, stationary fluid conflguration in the Kerr metric with 
flxed values of BH mass M and BH spin a parameters^ . The disk is symmetric under 
reflection in the equatorial plane. Matter is supplied at a steady rate, M, through a 
boundary "at infinity" and angular momentum is removed through the same boundary, 
whereas zero torque is assumed at the BH horizon. We include no self-irradiation of the 
disk, and we neglect the magnetic pressure, as well as the angular momentum carried 
away by radiation. Dissipation and angular momentum transport are given by the a 
prescription (Section 2.2.1), with a constant value of a in radius. We assume that the 
dissipation rate is proportional to the total pressure, p = + Prad- 

We take G = c = 1 and put r for the cylindrical radius. We also define the following 
expressions involving the BH spin: 

A = M\rl-2r, + al), (4.3) 

A = M\rt + rlal + 2nal), 

C = 1 - Sr^i + 2a,r;-^/2^ 

V = l-2r;' + 2alr;\ 

n = l-4a,r;3/2 + 3a2r;^ 

with a* = a/M and r* = r/M. The radial gas velocity, V, as measured by the observer 
co-rotating with the fiuid at a fixed value of r, is given by the relation (Abramowicz 
ET AL., 1996), 

V/VT^ = u^gli' = (4.4) 

We adopt an equation of state corresponding to the choice pgas = kpT/{fimp), and 
Prad = clT'^/S, with k the Boltzmann constant, rrip the proton mass, and a the radiation 
constant (no confusion with the spin parameter may arise). The mean molecular weight 
is taken to be /i = 0.62 corresponding to X = 0.7 and Y = 0.3. 

Only for the purpose of integrating the vertical structure, we assume the polytropic 
equation of state with the polytropic index N: p = Kp^^^^^ (N = 3 and K = const 
correspond to a mixture of gas and radiation with constant pressure ratio) . The vertical 



is the angular momentum per unit mass, a = J/M. while a» = a/AI = J/Al"^ is called the 
dimensionless spin parameter. 



4.1 Equations 



43 



integration of the hydrostatic equihbrium formula (Eq. 3.37), gives 

IP 

From the previous equation and the polytropic equation of state, it follows 



P = Po(l-^) . (4.5) 



T = To{l-i,\, (4.6) 



m J 

where po and Tq denote density and temperature at the equatorial plane, respectively. 
Equations 

In this section we present equations describing relativistic, stationary slim disks. They 
were derived originally by Lasota (1994) and improved by several other authors in- 
cluding Abramowicz et al. (1996), Gammie & Popham (1998) and Kato et al. 
(2008). More general derivation for the non-stationary case is given in Appendix B. 

(i) The mass conservation: 

M = -^.^y.A^^ — ^ (4.7) 

where S = p dz is disk surface density. 

(ii) The radial momentum conservation: 

V dV A IdP 

(4.8) 



where 



V"^ dr r S dr 



MA 



and fl = u'^/ M* is the angular velocity with respect to the stationary observer, Cl = fl — 
is the angular velocity with respect to the inertial observer, = ±M^/^/(r'^/^ ± 
aM^/^) are the angular frequencies of the co-rotating and counter- rotating Keplerian 
orbits, R = A/(PA^^'^) and P = pdz is vertically integrated pressure, 
(iii) The angular momentum conservation: 

— {C - An) = -^aP (4.10) 

zvr r 

where C = Urp, £in is the angular momentum at the disk inner edge, and 



"'=1— +— I f^") 



is the Lorentz factor. 

(iv) The vertical equilibrium: 
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Integration of Eq. 3.37, taking into account Eq. 4.5, leads to 



h^nl = {2N + 3)^, (4.12) 



where 



is the vertical epicyclic frequency and h corresponds to the disk thickness. Herein, we 
always assume = 3. 

(v) The energy conservation: 
The advective cooling is defined, in terms of the vertically integrated quantities, as 
(Kato ET AL., 2008) 



Q-^- = _ fru'{E + P))-u'—- u'^dz, (4.14) 
r dr dr oz 



where E is the height-integrated energy density. 



h 



-h 



5/3-1 



E= / 77?^ + 3Prad dz. (4.15) 



Using the mass conservation (Eq. 4.7), hydrostatic equilibrium (Eq. 3.37), and the 
following derivation. 



[ u''^dz = VL\ [ u'^pzdz 
J-h dz 

\nlf}lirpu^).^dz 
^ n2 d / 1 2 



fil— - / pz'dz 



2'Kr dr V S 
M ^2 dr/4 



2'Kr dr 



Q may be expressed as 



^,dv M / P dlnP ^ ^ ^P dlnS ^ 

27rr2 \Sdlnr Sdlnr 

Pdlnr/3 2 dlnr/4\ 

+ ^3^-7-^ + fi^r/4— . (4.16) 

L dmr dmr / 

Formulae for coefficients 773 and 774 are given below. 

The amount of heat advected Q'^^'^ is equal to the difference between viscous heating 
and radiative cooling: 

nadv ^vis ^rad ^ p^^^ S^dT^ 
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where Tc is the temperature at the equatorial plane. The opacity coefficient k is 
calculated using the Kramer's approximation (Eq. 3.20). 

(vi) The equation of state (in terms of vertically integrated quantities) 

P = + IviaT^- (4.18) 

/imp 3 

The radial derivative of the height-integrated pressure takes the form 

^ = (4-3/3)^ + /?^ + 
dr dr dr 

with (3 = ri2{k/ fimp)T,Tc/ P . 

The coefficients r]i to 774 are given by, and equal to (assuming Eqs. 4.5 and 4.6), 



1 '■^ 



7-14 

-'o ^0 



ciz = hh, (4.20) 



/" pTdz = In+i/In, (4.21) 

JO 



^1 
^2 



r/3 = - --^ETc + 2J4aT^/i , (4.22) 

P - 1 H In J 

1 

r/4 = - / pz^ dz = J^h^, (4.23) 



1 / 1 Jat+i 4 



where^ 



'AT 



a/tF r(l + A^) N&n 



2 r(3/2 + iV) (2Ar + l)!^ 



(4.24) 



_i r(3/2 + 7V) _ 1 

- 4r(5/2 + iV) - 6T4iV- ^^-^^^ 

4.2 Numerical methods 

The equations given above form a two-dimensional set of ordinary differential equations 
with a critical (sonic) point. One may choose e.g., V and Tc as the dependent variables. 
By a series of algebraic manipulations of Eqs. 4.7, 4.8, 4.10, 4.16, and 4.19, we obtain 
the following set of ordinary differential equations for V{r) and Tc(r)^: 

1 ^'""^ (4,26) 



r_128. T /T — & ■ 7 — 1 
^4 - 3T5' ^4/-'3 - 9, J3 - 18- 



^For clarity and compatibility with Chapter 6, we leave the radial derivatives of rji to 7^4 in Eqs. 4.27 
and 4.28. They are reduced to the derivatives of V and Tc, but the ultimate formulae are long and 
unreadable. 
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Mr^g P d In r d In r 

dln?73 fi^S ?74 dln?74 
dlnr P ?73 dlnr ' 



(4.27) 



with A/", V and Fi given by 



. 27rr2 ,^ P /r(r-M)~- 
+ ilj^V^i-^, (4.28) 



dlnr / ?73 dlnr 

V = ^'-fi|, (4.29) 

fi = 1 + (4.30) 

^3 

To obtain the regular solution one has to solve these equations together with outer 
boundary conditions given at some large radius rout and the following regularity condi- 
tions at the sonic radius rgon^ 

ATI-, =Vl-r =0 (4.31) 

' ' — ' son ' ' — ' son ^ ' 

which provides additional condition for calculating Ci^ which is the eigenvalue of the 
problem. The location of the sonic point is not known a priori. A mathematical problem 
defined in such a way can be classified as a two point boundary value problem which we 
solve applying the relaxation technique (PRESS, 2002). To start the relaxation process 
one has to provide a trial solution. The convergence strictly depends on the quality 
of such initial guess. In this work we apply the following method of searching for the 
proper initial conditions for relaxation. 

The angular momentum at the inner edge of the disk, £in, is the eigenvalue of 
the problem and must be chosen properly to satisfy the regularity conditions given by 
Eq. 4.31. The value of determines the shape of a solution. It turns out that the 
topology of slim disk solutions with higher and lower than the proper value £in,o 
is different. The former branch of solutions terminates as soon as the denominator V 
vanishes (the regularity condition is not satisfied). The latter does not cross the sonic 
point at all {T> is always positive). The self-consistent solution is expected to be the 
common limit of these two branches. A trial solution can be achieved in a few iteration 
steps. We start integration at r > lOOOM assuming the NoviKOV & Thorne (1973) 
boundary conditions and we use the implicit Runge-Kutta method of the 4th order. 

Once the trial solution is found one can start relaxation process between rout and the 
estimated position of the sonic radius rgon- The standard approach has to be modified 
due to the fact that we expect singularity at the inner boundary. Thus, we treat 
the problem as a free boundary problem and introduce one more variable describing 
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r [M] 

Figure 4.1: Temperature profiles in the vicinity of the sonic point for a few iterations 
leading to the trial model used as initial condition in the relaxation procedure (for 
details see Section 4.2). Solutions with too high value of terminate before reaching 
the proper sonic point while solutions with too low C\n go through the sonic point radius 
but follow an improper branch. The relaxed solution and the location of the sonic point 
are also presented. Models calculated for a non-spinning BH with M = O.lMEdd- 
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the position of the critical point (PRESS, 2002). In this work we usually use 100 mesh 
points spaced logarithmically in radius between the sonic point and the outer boundary. 

A relaxed solution is obtained in a few iteration steps and is then used as the 
initial condition for relaxation when looking for the solution of a problem with one 
parameter (e.g., mass accretion rate M or BH spin a) slightly modified. For the new 
system parameters we look for the outer boundary conditions at rout by integrating the 
equations from l.lrout (assuming the same outer boundary conditions as before) until 
we reach Tout- In such a way a full spectrum of system parameters can be achieved 
effectively. 

Once the solution outside the sonic point is found we numerically estimate the radial 
derivatives of V and Tc at the sonic point using values given at r > rson- Taking them 
into account we make a small step from the innermost mesh point (which corresponds 
to the sonic point) inward. Then we start integrating using standard Runge-Kutta 
method until we get close enough to the horizon. 



4.3 Solutions 
Definitions 

We introduce the following critical (Eddington) luminosity, 



= 1^.1. 25x10-^ erg.-', (4.32) 



and mass accretion rate^, 

Mem = 16 X ^ = 16 x = 2.23 xlO^« t^S^ , 4.33 



where /tes = 0.20(1 + X) cm^/g is the electron scattering opacity and here we take X = 1 
being the hydrogen mass fraction. We also define the efficiency of accretion as 

in ^ ^1 -^Edd f A 

?7 = 16 X —. = — : . (4.34) 

Mc2 M/MEdd 

The factor 16 in Eqs. 4.33 and 4.34 resembles the fact that the efficiency of accretion in 
the Pseudo-Newtonian potential for a thin, Shakura-Sunyaev disk, is 1 — n^ in = 1/16 
(compare Section 2.3). In such a way we may consider MEdd as an accretion rate for 
which, in case of a non-rotating BH only, the disk luminosity is close to -^Edd- We also 
introduce a dimensionless accretion rate rh, 

M 

m = -. . 4.35 

MEdd 



*Many authors use a different definition, M^^^^ — L^dd/c^- 
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r [M] 

Figure 4.2: Flux profiles for different mass accretion rates in case of a non-rotating BH 
and two values of a: 0.01 (black solid), 0.1 (red dashed lines). For each value five lines 
for the following mass accretion rates: 0.01, 0.1, 1.0, 2.0 and lO.OMEdd are presented. 
BH mass is lOM©. 

Disk appearance 

In this paragraph we discuss most important properties of slim disk solutions for the 
purpose of spectra modeling. 

In Fig. 4.2 we present profiles of the emitted flux for a Schwarzschild BH and a 
range of accretion rates. Profiles for two values of the a parameter are shown: 0.01 
(solid black) and 0.1 (red dashed lines). In the limit of low accretion rates the flux 
profiles coincide with the Novikov & Thorne solution - no flux emerges from inside 
the marginally stable orbit. The higher the accretion rate, the higher departure from 
the standard profile as the advection of heat becomes more important. For accretion 
rates exceeding MEdd the disk emission is significantly shifted inward. The emission 
profile in the inner disk regions for super-Eddington accretion rates has a different 
slope (F oc r~°'^) instead of oc r~°'^^ for radiatively efficient disks. Such behavior is 
valid for all values of the viscosity parameter a. However, the higher the value of a, 
the earlier advection affects the disk emission. 

Fig. 4.3 presents profiles of emission from a disk with m = 0.1 for a few values of 
BH spin. They all coincide at large radii, where the impact of BH rotation is negligible. 
The larger BH spin, the closer to the horizon the marginally stable orbit is, and, as a 
result, the disk emission extends much closer to the horizon, well below r = 6M. It is 
clear, that the luminosity of a disk increases with BH spin. This increase of luminosity 
corresponds exactly to the increase of the efficiency of accretion (Eq. 2.30) for low 
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Figure 4.3: Flux profiles for m = 0.1 and different values of BH spin. Results for two 
values of a: 0.01 (black solid), 0.1 (red dashed lines) are presented. BH mass is IOMq. 

accretion rates. 

Profiles of the radial velocity (K) and surface density (S) are shown in the top and 
bottom panels of Fig. 4.4, respectively. The black lines present solutions for a = 0.01, 
while the red dashed lines for a = 0.1. The radial velocity for large radii coincides with 
values given by NOVIKOV & Thorne (1973) which are assumed as the outer boundary 
conditions, and approaches the speed of light when getting close to the horizon. At a 
given radius the radial velocity increases with mass accretion rate. The radial velocity is 
directly related to the surface density through the mass conservation equation (Eq. 4.7). 
For large radii, where the disk is gas pressure dominated, the surface density increases 
with accretion rate (this region corresponds to the lower branch of solutions on the 
(T,E) plane, (see e.g., AbramowiCZ ET AL. (1988)) while for moderate radii the 
relation is opposite (the middle branch). For the highest accretion rates (rh > 10) 
the solutions enter the upper, advection dominated branch. As the surface density is 
expected to be inversely proportional to a, its values for a = 0.1 are roughly ten times 
smaller than for a = 0.01. 

It has been proven (S4DOWSKI ET AL., 2009) that the disk photosphere location 
should be taken into account when performing ray-tracing of photons emitted from an 
accretion disk. It becomes of major importance when the accretion rates are high and 
disks are no longer geometrically thin. In Fig. 4.5 we plot profiles of the photosphere 
(disk height, cosQh = -f^phot /''"spherical) for different accretion rates and a non-rotating 
BH. It is clearly visible that for accretion rates m > 0.1 the inner regions are thicker due 
to the radiation pressure. For M > LOMedd the highest value of cos Qh ratio exceeds 
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Figure 4.4: Profiles of the radial velocity (top panel) and surface density (bottom panel) 
for slim disk solutions with different accretion rates in case of a non-rotating BH with 
-^BH = 1OM0. Solid black lines present solutions for a = 0.01 while red dashed lines 
for a = 0.1. 
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Figure 4.5: Profiles of tlie disk pliotospliere in terms of cos©// = i^phot /'"sph for a 
number of accretion rates and two values of a = 0.01 and 0.1. Mbh = IOMq and 
a* = 0. 



0.3. In Fig. 7.1 similar profiles are plotted using the same scale on both axes. 

In Fig. 4.7 we present profiles of the angular momentum for two representative 
values of a (0.01 and 0.1 on left and right panel, respectively) and a few accretion 
rates. For low accretion rates the angular momentum follows the Keplerian profile 
(down to the marginally stable orbit) for both cases. For higher accretion rates it 
departs from the Keplerian flow. For the lower value of a the flow is super-Keplerian 
for moderate radii. The higher mass accretion rate the more significant deviation from 
the Keplerian profile. For the higher a (0.1) the behavior is quantitatively different 
as the angular momentum at the horizon (jCin) falls below the Keplerian value at the 
marginally stable orbit. The higher the accretion rate the smaller the value. The impact 
of a on the slim disk solution topology has been extensively studied in Abramowicz 
ET AL. (2010). Figure 4.8 shows how £in depends on both the accretion rate and the 
a viscosity parameter. 

Accretion disks are expected to be gas-pressure dominated in the limit of low accre- 
tion rates (Shakura & SUNYAEV, 1973). However, the temperature, and correspond- 
ing radiation pressure, increases with accretion rate. As a result, for accretion rates 
higher than a particular value (e.g., O.lMEdd for Ci* = 0), there is a radiation-pressure 
dominated region, the extent of which increases with further grow of accretion rate. 
This behavior is presented in Fig. 4.9 which shows the gas pressure to total pressure 
ratio (/3) for different accretion rates and two representative values of a. For m = 0.1, 
the region between 8M and 50M is dominated by radiation pressure. The pressure 
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Figure 4.6: Disk height profiles of slim disks with three accretion rates: 0.1 (blue), 1.0 (orange) and lO.OMEdd (green line), and 
a = 0.1. The black circle denotes the central BH. 
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Figure 4.7: Profiles of tlie disk angular momentum (u^) for a = 0.01 (top) and a = 0.1 
(bottom panel) for different accretion rates. The spin of the BH is a* =0. 



4.3 Solutions 



55 




Figure 4.8: Dependence of the angular momentum at the horizon on accretion rate for 
solutions with different values of a for = 0. 



ratio drops below 0.4 which is the critical value for triggering the thermal instability 
according to the standard theory (for detailed discussion see Section 5.1). For the Ed- 
dington accretion rate (m = 1.0) the unstable region extends between 5M and 200M. 
For super-Eddington accretion rates the disk is radiation-pressure supported down to 
BH horizon. 

Standard thin accretion disks are radiatively efficient — viscous heating is balanced 
by radiative cooling. This is not the case for accretion rates close to, and higher than 
the Eddington one. Such disks are thick and are characterized by high radial velocity 
of gas. As a result, additional mechanism of cooling triggers in: advection. In Fig. 4.10 
we plot the advection coefficient f'^'^^ defined as the ratio of the advective to radiative 
cooling rates. Values close to zero denote radiatively efficient disk regions, higher than 
unity — regions where advective cooling dominates over radiative cooling, negative — 
regions where heat advected from outer disk annuli is released (advective heating). Disk 
solutions for the lowest accretion rates are in fact radiatively efficient (/^*^^ ~ 0). For 
accretion rates close to MEdd the region r > lOM is significantly affected by advection, 
e.g., the advective cooling equals up to 30% of the radiative cooling for m = 2.0. 
The region inside r = lOM is characterized by negative advection coefficient. Radiative 
cooling balances both viscous and advective heating. The heat advected into this region 
is the only heating mechanism inside ISCO {f^'^^ ~ ~1); causing the emission from 
inside this radius visible in Fig. 4.2. For strongly super-Eddington accretion rates (e.g., 
rh = 10) the advective cooling rate overwhelms the rate of radiative cooling outside 
ISCO, while advection heating still dominates close to BH horizon. 

The advection of heat in slim disks significantly changes disk properties for lumi- 
nosities L > IvEdd- This fact is clearly visible in Fig. 4.11. In the top panel we plot the 
relation between disk luminosity and mass accretion rate for three values of BH spin. 
For sub-Eddington regime the luminosities are proportional to the accretion rates as the 
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Figure 4.9: Gas pressure to total pressure ratio for slim disk solutions with m = 0.01, 
0.1, 1.0 and 10.0. Profiles for a = 0.01 and 0.1 are presented with solid black and 
dashed red lines, respectively. 

efficiency of accretion (77 = 16 L/Mc^, plotted in the bottom panel) remains constant. 
For higher accretion rates, however, the efficiency drops down as some amount of heat 
generated by viscosity is captured in the flow and advected onto the BH. The higher 
the accretion rate, the lower the efficiency of accretion. For = 0.9 it decreases by an 
order of magnitude for m ^ 10. 

In Fig. 4.12 we present a schematic picture of a slim accretion disk with super- 
Eddington accretion rate. 

4.4 Inner edge(s) 

Accretion flows onto BHs must change character before matter crosses the event horizon, 
for two reasons. First, matter must cross the black- hole surface at the speed of light 
as measured by a local inertial observer (GOURGOULHON & Jaramillo, 2006), so 
that if the flow is subsonic far away from the black-hole (in practice it is always the 
case) it will have to cross the sound barrier (well) before reaching the horizon. This is 
the property of all realistic flows independent of their angular momentum. This sonic 
surface can be considered as the inner edge of the accretion flow. 

The second reason is related to angular momentum. Far from the hole many (most 
probably most) rotating accretion flows adapt the Keplerian angular momentum profile. 
Because of the existence of the innermost stable circular orbit (ISCO), these flows must 
stop to be Keplerian there. At high accretion rates when pressure gradients become 
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Figure 4.10: The advection factor profiles for rh = 0.01, 2.0 and 10.0. Profiles for 
a = 0.01 and 0.1 are presented with solid black and dashed red lines, respectively. The 
fraction f^'^^/{l + J^*^^) of heat generated by viscosity is accumulated in the flow. In 
regions with f^'^^ < the advected heat is released. 
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Figure 4.11: Top panel: Luminosity vs accretion rate for three values of BH spin 
(a* = 0.0, 0.9, 0.999) and two values of a = 0.01 (black) and 0.1 (red line). Bottom 
panel: efficiency of accretion t] = 16 L/Mc^. The luminosity L was integrated locally, 
not taking into account the relativistic g-factor (see Section 4.5.2), and therefore only 
roughly corresponds to fj defined in Eq. 2.30. 




Figure 4.12: A scheme of an accretion disk at super-Eddington accretion rate: Gas is accreted onto a BH through a disk that 
significantly extends above the equatorial plane (geometrically slim). Viscosity results from turbulent motion. The profile of 
rotation departs from the Keplerian one. Disk extends down to BH horizon. Viscous heating is balanced by radiative and 
advective cooling mechanisms. The disk emission does not terminate at ISCO — significant emission comes from the vicinity of 
the BH, resulting from heat advected from outer regions. 
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important, the flow may extend below the ISCO but the presence of the innermost 
bound circular orbit (IBCO) defines another limit to a circular flow (the absolute limit 
being given by the circular photon orbit, the CPO). These critical circular orbits provide 
another possible definition of the inner edge of the flow, in this case of an accretion 
disk. 

The question is: what is the relation between the accretion flow edges? In the case 
of geometrically thin disks the sonic and Keplerian edges coincide and one can define 
the ISCO as the inner edge of these disks. Paczynski (2000) showed rigorously that 
independent of e.g., viscosity mechanism and the presence of magnetic fields, the ISCO 
is the universal inner disk's edge for not too-high viscosities. The case of thin disks is 
therefore settled^. 

However, this is not the case for non-thin accretion disks, i.e., of medium and high 
luminosities. The problem of defining the inner edge of an accretion disk is not just 
a formal exercise. Afshordi & PACZYNSKI (2003) explored several reasons which 
made discussing the precise location of inner edge r = r^^ of the BH accretion disks an 
interesting and important issue. One of them was. 

Theory of accretion disks is several decades old. With time ever more sophisticated 
and more diverse models of accretion onto BHs have been introduced. However, when 
it comes to modeling disk spectra, conventional steady state, geometrically thin-disk 
models are still used, adopting the classical "no torque" inner boundary condition at 
the marginally stable orbit. 

Krolik & Hawley (2002) proposed several "empirical" definitions of the inner 
edge, each serving a different practical purpose (see also the follow-up investigation by 
Beckwith ET AL. (2002)). We add to these a few more definitions. The hst of the 
inner edges considered in this section consists of^, 

1. The potential spout edge r^^ = '"pot, where the effective potential forms a self- 
crossing Roche lobe, and accretion is governed by the Roche lobe overflow. 

2. The sonic edge Tin = rson, where the transition from subsonic to transonic accre- 
tion occurs. Hydrodynamical disturbances do not propagate upstream a super- 
sonic flow, and therefore the subsonic part of the flow is "causally" disconnected 
from the supersonic part. 

3. The variability edge = r^ar, the smallest radius where orbital motion of coherent 
spots may produce quasi-periodic variability. 

^Penna et AL. (2010) studied the effects of magnetic fields on thin accretion disk (the disk 
thickness h/r < 0.07. They found that to within a few percent the magnetized disks are consistent 
with the NoviKOV & Thorne (1973) model, in which the inner edge coincides with the ISCO. 

^Krolik & Hawley defined the inner edges no. 4, 5 and 6 and in addition a seventh edge 7, the 
turbulence edge, where flux-freezing becomes more important than turbulence in determining the mag- 
netic field structure. Magnetic fields are not considered for slim accretion disks, and we therefore do 
not discuss 7. 
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4. The stress edge Hn = rgtr, the outermost radius where the Reynolds stress is small, 
and plunging matter has no dynamical contact with the outer accretion flow. 

5. The radiation edge rin = Trad, the innermost place from which significant lumi- 
nosity emerges. 

6. The reflection edge r^^ = Tref , the smallest radius capable of producing significant 
fluorescent iron line. 

In the next six subsections, we discuss these six edges one by one basing on the 
solutions of the relativistic slim disk model presented in the previous section. This 
discussion is similar to the study of Abramowicz ET AL. (2010) but is based on a slim 
disk model with more consistent (i.e., polytropic) treatment of the vertical structure. 
Qualitative behavior is exactly the same. However, the impact of advection sets on 
at a higher (0.6 versus O.SMEdd) accretion rate. We also point out here that locations 
of some of these inner edges are sensitive to the assumed value of the a parameter. 
Solutions presented in the previous Secion are based on the a = const assumption. 
This assumption, however, may not be satisfied inside the marginally stable orbit, as 
some of the most recent GRMHD simulations show (e.g., Penna ET AL. (2010)). 
Results presented in this Section should be therefore considered qualitative only. 

4.4.1 Potential spout edge 

The idea of the "relativistic Roche lobe overflow" governing accretion close to the BH 
was first explained by Paczyhski (see KOZLOWSKI ET AL. (1978)). It was later explored 
in detail by many authors analytically (e.g., Abramowicz (1981, 1985)) and by large- 
scale hydrodynamical simulations (e.g., IGUMENSHCHEV & BELOBORODOV (1997)). It 
became a standard concept in BH accretion theory. Figure 4.13 schematically illustrates 
the Roche lobe overflow mechanism. The leftmost panel presents a demonstrative 
proflle of disk angular momentum, which reaches the Keplerian value at the radius 
corresponding to the self-crossing of the equipotential surfaces presented in the middle 
panel. To flow through this "cusp", matter must have potential energy higher than the 
value of the potential at this point, i.e., the "potential barrier" is crossed only when the 
matter overflows its Roche lobe. 

The potential difference between the horizon and the spout is inflnite, and therefore 
no external force can prevent the matter located there from plunging into the BH. At 
radii greater that rpot, the potential barrier at r = rpot retains the matter inside this 
radius. We note, that because the dynamical equilibrium is given (approximately) by 
ViWcff = VjP/ p, with p being the density, one may also say that it is the pressure 
gradient (the pressure stress) that holds the matter within rpot- 

The speciflc angular momentum in the Novikov-Thorne model is assumed to be Ke- 
plerian. Slim disk models do not a priori assume an angular momentum distribution, 
but self-consistently calculate it from the relevant equations of hydrodynamics (see Sec- 
tion 4.1). These calculations indicate that the type of angular momentum distribution 
depends on whether the accretion rate and viscosity constrain the flow to be either 
disk-like or Bondi-like type. 
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Figure 4.13: An illustrative visualization of the Roche lobe overflow. The leftmost panel schematically presents disk an- 
gular momentum profile and its relation to the Keplerian distribution. The middle panel shows the equipotential sur- 
faces. The dotted region denotes the volume filled with accreting fluid. The rightmost panel presents the potential bar- 
rier at the equatorial plane (z = 0) and the potential of the fluid {Ws) overflowing the barrier. The figure is taken from 
http : / /www . scholarpedia . org/ article/Accretion_discs. 
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In the Bondi-type accretion flows, the angular momentum is everywhere sub-Keplerian, 
C < Ck- These flows are typical of high viscosities and high accretion rates, as the 
case of a = 0.1 and rh = 100 shown in the bottom panel of Fig. 4.7. This is the only 
Bondi-like flow in this flgure. 

In the disk-like accretion flows, the angular momentum of the matter in the disk 
is sub-Keplerian everywhere, except the strong-gravity region rpot < r < Vccn, where 
the flow is super- Keplerian, C > Ck- The radius Vcen > fisco corresponds to the 
ring of the maximal pressure in the accretion disk. This is also the minimum of the 
effective potential. The radius rpot < ''"isco delineates a saddle point for both pressure 
and effective potential; this is also the location of the "potential spout inner edge". 



We note that in the classic solutions for spherically accretion flows found by BONDI 
(1952) the viscosity is unimportant and the sonic point is saddle, while in the "Bondi- 
like" flows discussed here, angular momentum transport by viscosity is essentially im- 
portant and the sonic point is usually nodal. Therefore, one should keep in mind that 
the difference between these types of accretion flows is also determined by the relative 
importance of pressure and viscosity. For this reason, a different terminology is often 
used. Instead of "disk-like", one uses the term "pressure-driven", and instead of "Bondi- 
like", one uses "viscosity-driven" (see e.g., Matsumoto et AL. (1984); Kato et AL. 



From the above discussion, it is clear that the location of this particular inner edge 
Tpot is formally given as the smaller of the two roots, r± = (r_|_, r_), of the equation 



The larger root corresponds to rccn- Obviously, Eq. 4.36 has always a solution for the 
disk-like flows, and never for the Bondi-like flows. 

The location of the potential spout inner edge rpot is shown in Fig. 4.14 for a = 0.01. 
We note that for low accretion rates, (e.g., rh < 0.6 for a* = 0), the location of 
the potential spout inner edge coincides with ISCO. At m ^ 0.6, the location of the 
potential spout jumps to a new position, which is close to the radius of the innermost 
bound circular orbit, Tibco = 4M. This behavior has long been recognized flrst by 
KOZLOWSKI ET AL. (1978) for Polish doughnuts, and then by Abramowicz ET AL. 
(1988) for slim disks. 

4.4.2 Sonic edge 

As it has already been discussed, by a series of algebraic manipulations, one may reduce 
the slim disk equations into a set of two ordinary differential equations for two dependent 
variables, e.g., the radial velocity V and the central temperature Tc, 



^in '"pot' 



(2008)). 



C{r) - CK{r) = 0. 



(4.36) 



dlnr 

dlnTc 
dlnr 



'D{r,V,Tc) ' 
Af2{r,V,Tc) 
'D{r,V,Tc) ■ 



(4.37) 
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Figure 4.14: Location of the potential spout inner edge rpot for viscosity a = 0.01 and 
a* = 0. Solid lines show the exact location of rpot given by Eq. (4.36). 

For a non-singular physical solution, the numerators A/i and A/2 must vanish at the 
same radius as the denominator V. The denominator vanishes at the sonic edge (or 
sonic radius) where the Mach number is close to unity, i.e. 

V{r,V,Tc)\r=r.^^ = 0. (4.38) 

For low mass accretion rates, lower than about O.GMEdd in the case of = 0, the sonic 
edge Tson is close to ISCO, independently of the viscosity a, as Fig. 4.15 shows. At 
about O.GMEdd; a qualitative change occurs, resembling a "phase transition" from the 
Shakura-Sunyaev behavior to a very different slim disk behavior. 

For higher accretion rates, the location of the sonic point significantly departs from 
ISCO. For low values of a, the sonic point moves closer to the horizon down to ~ 4M for 
a = 0.001. For a > 0.3, the sonic point moves outward with increasing accretion rate 
reaching values as high as 7.5M for a = 0.5 and lOOMEdd- This effect was first noticed 
for low accretion rates by MuCHOTRZEB-CzERNY (1986) and later investigated for a 
wide range of accretion rates by Abramowicz ET AL. (1988), who explained it in 
terms of the disk-Bondi dichotomy. The dependence of the sonic point location on the 
accretion rate in the near-Eddington regime is more complicated and is related to, for 
this range of accretion rates, the transition from the radiatively efficient disk to the 
slim disk occurring close to the sonic radius. It is also sensitive to the coefficients of the 
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Figure 4.15: Location of the sonic point as a function of the accretion rate for different 
values of a, for a non-rotating black hole, a* = 0. 



vertical integration introduced in Eqs. 4.20 — 4.23, as comparison with similar plots in 
S4DOWSKI (2009) and Abramowicz ET AL. (2010) shows. 

The topology of the sonic point is important, because physically acceptable solutions 
must be of the saddle or nodal type, the spiral type being forbidden. The topology may 
be classified by the eigenvalues Ai, A2, A3 of the Jacobi matrix 

r dV dV dV_ 

dr dV dTc 

rr _ aM dhh dMl 

^ — dr dV dTc 

dAf2 8X2 dN2 
L dr aV dTc 

Because det(jr) = 0, only two eigenvalues Ai, A2 are non-zero, and the quadratic char- 
acteristic equation that determines them takes the form, 

2 A^ - 2 A ii{J) - {ii{J^) - tT\j)) = 0. (4.40) 

The nodal type is given by A1A2 > and the saddle type by A1A2 < 0. For the 
lowest values of a, only the saddle-type solutions exist. For moderate values of a 
(0.1 < a < 0.4), the topological type of the sonic point changes at least once with 
increasing accretion rate. For the highest a, solutions have only nodal-type critical 
points. 



(4.39) 



4.4.3 Variability edge 

Axially symmetric and stationary states of slim accretion disks are, obviously, theoret- 
ical idealizations. Real disks are non-axial and non-steady. In particular, one expects 
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transient coherent features at accretion disk surfaces — clumps, flares, and vortices. 
The orbital motion of these features could quasi-periodically modulate the observed 
flux of radiation, mostly by means of the Doppler effect and the relativistic beaming. 
We define 11 to be the "averaged" variability period, and All a change in the period 
during one period caused by the radial motion of a spot. The variability quality factor 
Q may be estimated by, 

1 An Af] 1 dl] , 1 dfi , 

oTJT^^ = ^^7^.-r--7. (4.41) 



g n n vtdr n"^ dr u 

where u'' /u^ = dr/dt and u"^' and are contravariant components of the four velocity. 
The period is given in terms of the orbital angular velocity by 11 = 27r/r2. Using the 
relations (see Eqs. 4.3 for the explanation of the notation used) 

V 

u 



VI - r 



= ^ = ^ , (4.42) 



rVA 



{l-V^){l-{V'i>Y) 



where V is the radial velocity measured by an observer co-rotating with the fluid, one 
obtains 

1 dlogl] 

Q = ir -r^ 17 / (a*,r), (4.43) 

Itx d log r V 

where 

.5„4^-l/2 

- y\ - y\ a^ya^^ L) - y\ 

Y4> 



r(a„r) ^ -—={i-X-X^al{al + l)-X'a':)~'^\ (4.44) 



= M, 

and X = 2M/r. From their definitions, it is clear that AA > outside the BH horizon. 
We note that in the Newtonian limit it is X -C 1 and one has /*(«*, r) = 1. In this 
limit, V and V'l' are the radial and azimuthal components of velocity, and Eq. 4.43 
takes its obvious Newtonian form. 

The behavior of the quality factor Q is shown in Fig. 4.17. Profiles for four accretion 
rates are drawn. As Fig. 4.16 shows, the lower accretion rate, the smaller radial velocity 
component, hence the quality factor Q in general increases with decreasing accretion 
rate. For the lowest values of rh, a rapid drop is visible at ISCO corresponding to the 
change in the nature of the flow (gas enters the free-fall region below ISCO). For higher 
accretion rates, this behavior is suppressed as the trajectories become wide open spirals 
well outside ISCO. 

We note that our definition given in Eq. 4.41 of the quality factor Q, essentially 
agrees with a practical definition of the variability quality factor Qq defined by observers 




Figure 4.16: The fluid flow trajectories in slim accretion disks shown by thin solid lines for different accretion rates. Locations 
of Tpot and the location of black hole horizon are shown by thick orange solid and broken lines, respectively. For low accretion 
rates, the pattern of trajectories consists of very tight spirals (almost circles) for r > Tpot ~ ''"isco and very wide spirals (almost 
a radial fall) for r < Tpot- In this case, there is a sharp transition from almost circular motion to almost radial free-fall that 
clearly defines the variability edge as r^ar = '"isco- For higher accretion rates, the fluid trajectories are wide open spirals in the 
whole inner region of the flow and the variability edge makes no sense. 
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Figure 4.17: The quality factor Q profiles for different accretion rates and a = 0.01. 
Triangles show rpot for each rate. The vertical dashed line denotes the location of ISCO. 

with the help of the observationally constructed Fourier variability power spectra, /(z^). 
Here I{u) is the observed variability power (i.e. the square of the observed amplitude) at 
a particular observed variability frequency u. Any observed quasi periodic variability 
of the frequency ~ uq can clearly be seen in the power spectrum as a local peak in 
J(i^), centered on a particular frequency uq. The half-width Au of the peak defines the 
variability quality factor by Qq = u/Auq. 

Quasi-periodic variability of kHz frequencies, called kHz QPO, is observed from 
several low-mass neutron star and BH binaries. In a pioneering and important piece of 
research, BARRET ET AL. (2005) carefully measured the quality factor for a particular 
source in this class (4U 1608-52) and found that Qq ~ 200, i.e. that the kHz signals are 
very coherent. They argued that Qq ~ 200 cannot be due to kinematic effects in the 
orbital motion of hot spots, clumps, or other similar features located at the accretion 
disk surface, because these features are too quickly sheared out by the differential 
rotation of the disk (see also Bath et AL. (1974); Pringle (1981)). They also 
argued that although coherent vortices may survive much longer times at the disk 
surface (e.g., Abramowicz et AL. (1995)), if they participate in the inward radial 
motion, the observed variability Qq cannot be high. Our results shown in Fig. 4.17 
illustrate and strengthen this point. We also agree with the conclusion reached by 
Barret et AL. (2005) that the observational evidence against orbiting clumps as 
a possible explanation of the phenomenon of kHz QPO, seems to indicate that this 
phenomenon is probably caused by the accretion disk global oscillations^. For excellent 



''Barret et AL. (2005) also found how Qq varies in time for each of the two individual oscillations 
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reviews of the QPO oscillatory models, we refer to WAGONER (1999) and Kato (2001). 

Although clumps, hot-spots, vortices or magnetic flares cannot explain the coherent 
kHz QPOs with Qq ~ 200, they certainly are important in explaining the continuous, 
featureless Fourier variability power spectra (see e.g., Abramowicz et AL. (1991); 
SCHNITTMAN (2005); Pechacek, Karas & CzERNY (2008) and references quoted 
there). Our results shown in Fig. 4.17 indicate that: (i) at low accretion rates, a sharp 
high-frequency cut-off in /(i/) may be expected at about the ISCO frequency; (ii) at 
high accretion rates, there should be no cut-off in /(z/) at any frequency; and (iii) the 
logarithmic slope p = {dlnl /dlnu) should depend on rh. 

AAA Stress edge 

The Shakura-Sunyaev model assumes that there is no torque at the inner edge of the 
disk, which in this model coincides with ISCO. Slim disk model assumes that there is 
no torque at the horizon of the BH. It makes no assumption about the torque at the 
disk inner edge, but calculations prove that the torque is small there. 

The zero-torque at the horizon is consistent with the small torque at the inner 
edge of slim disks, as Fig. 4.18 shows. The Figure presents the relative importance 
of the torque T by comparing it with the "advective" flux of angular momentum Mj 
(J = Mj + T). For the viscosity parameter a smaller than about 0.01, the ratio 
X = T/Mj at both rpot and rgon is smaller than 0.01 even for highly super-Eddington 
accretion rates, and for low accretion rates the ratio is vanishingly small, x ~ 10~^. For 
high viscosity, a = 0.5, the ratio is small for low accretion rates, x < 10~^ and smaller 
than 0.1 for super-Eddington accretion rates (calculated at the sonic radius, as the disk 
enters the Bondi-like regime at these high accretion rates). For two highest viscosity 
models (a = 0.3 and 0.5) the x parameter calculated at rj„ exceeds 0.1 for accretion 
rates close to the disk-Bondi limiting value what may suggest the no-stress at the disk 
inner edge assumption is violated for such high viscosities. 

To define the stress inner edge r^tj., one has to specify the characteristic value of the 
torque parameter x- Profiles of rgtr for a few values of x and a = 0.01 are shown in 
Fig. 4.19. The stress edge for x — )■ is located at ISCO for low accretion rates. When 
the accretion rate exceeds ~ 0.6MEdd, the edge departs from ISCO and moves closer to 
BH approaching its horizon with increasing m. The behavior of rgtr profiles for higher 
(> 0.1) values of x is different - they move away from the BH as the angular momentum 
profiles become flatter with increasing accretion rates (compare Fig. 4.7). 

By pushing the MHD numerical simulations to their limits, Shafee ET AL. (2006) 
and Penna et AL. (2010) calculated a thin, H/r < 0.1, disk-like accretion flow, and 
demonstrated that its inner edge torque is in fact small. 

4.4.5 Radiation edge 

As discussed in the previous section, the torque at Vpot < tisco is small, but non-zero 
and therefore there is also orbital energy dissipation at radii smaller than ISCO. Thus, 

in the "twin-peak QPO". This places strong observational constraints on possible oscillatory models of 
the twin peak kHZ QPO; see also Boutelier et al. (2010). 
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Figure 4.18: Ratio of the angular momentum flux caused by torque to the flux caused 
by advection calculated at rpot (top) and Tgon (bottom panel) versus mass accretion rate 
for a number of values of a and a* = 0. The rpot profiles for high viscosities terminate 
when the disk enters the Bondi-like regime. 
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Figure 4.19: Profiles of r^tr defined as the radius with given value of the torque parameter 
X for a = 0.01. BH horizon and ISCO are also shown with dot-dashed and dashed lines, 
respectively. 

some radiation does originate in this region. Moreover, the heat advected with the 
flow may be radiated away in this region. As a result, the radiation edge Trad is not 
expected to coincide with ISCO. In Fig. 4.20, we present profiles of Trad defined as 
the radii limiting area emitting a given fraction of the disk total luminosity. For low 
accretion rates (< 0.6MEdd), the disk emission terminates close to ISCO as the classical 
models of accretion disks predict. The locations of the presented Trad are determined by 
the regular Novikov & Thorne flux radial profile. For higher accretion rates, the disk 
becomes advective and the maximum of the emission moves significantly inward. As a 
consequence of the increasing rate of advection (and resulting inward shift of Tpot), the 
efficiency of accretion drops down. 

We emphasize that the location of the radiation edge is not determined by the 
location of the stress edge (as some authors seem to believe), but by the significant 
advection flux bringing energy into the region well below ISCO. 



4.4.6 Reflection edge 

The iron Kq, fluorescent line is an observed characteristic feature of many sources with 
BH accretion disks (Miller, 2006; Remillard & McClintock, 2006). The in- 
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Figure 4.20: "Luminosity edges" defining the inner radii of tlie area emitting a given 
amount of tlie total disk radiation. Tlie lines are drawn for 95%, 99%, and 99.9% of 
the total emission. The dashed line shows the location of the potential spout inner edge 
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The gravitational suppression of the radiation has been taken into account. 



tensity and the shape of this line depends strongly on the physical conditions close to 
the inner edge. This has been discussed by many authors, including REYNOLDS & 
Fabian (2008) who gave three conditions for line formation: (i) the flow has to be 
Thomson-thick in the vertical direction; (ii) disk has to be irradiated by an external 
source of X-rays (hard X-ray irradiation plays a crucial role in the process of fluores- 
cence and changes the ionization degree of matter); and (iii) the ionization state should 
be sufficiently low (iron cannot be fully ionized). 

Nevertheless, since fluorescent iron-line emission depends on many aspects, such 
as the energy distribution of illuminating photons, temperature, ionization state, and 
density of the emitting matter as well as iron abundance, there is no obvious condition 
for the reflection edge (defined as the minimal radius where the refiection features 
originate). Additional computations of refiection models show that for some set of 
parameters iron fiuorescent line may arise even from Thomson-thin matter (DUMONT 
ET AL., 2002). In this paper, we assume that the refiection edge is connected to the 
condition 

TcS = A/Tabs(Tabs + ^s) > I- (4.45) 

However, one has to keep in mind that the effective optical depth at the iron line band 
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may be much larger than the above, frequency-averaged value. 

In Fig. 4.21, we show profiles of the effective optical depth Tcs in different regimes of 
accretion rates for a = 0.1 and a* = 0. Three characteristic types of their behavior are 
shown: sharp drop, maximum and monotonic at the top, middle, and bottom panels, 
respectively. The behavior of different values of a and is qualitatively similar (but 
not quantitatively as in general t^s increases with decreasing a). The top panel of the 
figure, corresponding to the lowest accretion rates, shows a sharp drop in near ISCO. 
The same behavior was noticed previously e.g., by REYNOLDS & FABIAN (2008). The 
drop may clearly define the inner reflection edge r^-d ~ ^'isco; limiting the radii where 
formation of the fluorescent iron line is prominent. The middle panel, corresponding 
to moderate accretion rates, shows a maximum in t^s near ISCO. The non-monotonic 
behavior is caused by the regions of moderate radii outside ISCO being both radiation 
pressure and scattering dominated. We note, that the top of the maximum of Teg 
remains close to ISCO for a range of accretion rates, but for accretion rates higher than 
O.GMEdd, it moves closer to the black hole with increasing rh as the disk emission profile 
changes due to advection. The bottommost panel corresponds to super-Eddington 
accretion rates. The profiles are monotonic in and define no characteristic inner 
reflection edge. Close to the BH, these flows are effectively optically thin reaching 
Teff = 1 at about a few tens of gravitational radii. 

When the effective optical depth of the flow becomes less then unity, opacities be- 
comes invalid. In these cases, full radiative transfer through accretion disk atmospheres 
has to be considered (e.g., DAVIS ET AL. (2005); RozANSKA & Madej (2008)). Nev- 
ertheless, our results allow us to estimate roughly how far from the black hole the iron 
line formation is most prominent, assuming that the disk is uniformly illuminated by 
an exterior X-ray source. For accretion rates lower than O.GMEdd, the reflection edge is 
located very close to ISCO and we may identify the shape of the iron line with the grav- 
itational and dynamical effects of the ISCO. In the case of higher but sub-Eddington 
accretion rates, the maximum of the effective optical depth is located inside the ISCO, 
which may possibly allow us to study extreme gravitational effects on the iron line 
profile. However, the assumption that the line is formed at the ISCO is no longer 
satisfied. The super-Eddington flows have smooth and monotonic profiles of effective 
optical depth. Therefore, the reflection edge cannot be uniquely defined and no relation 
between the shape of the fluorescent lines and ISCO exists. Finally, we note that these 
lines may be successfully modeled by clumpy absorbing material and have nothing to 
do with relativistic effects (see e.g., MiLLER ET AL. (2009) and references therein). 
The role of the ISCO in determining the shape of the Fe lines was also questioned in 
the past (based on different reasoning) by REYNOLDS & Begelman (1997), whose 
arguments were then refuted by YoUNG ET AL. (1998). 

4.4.7 Summary 

We have addressed the inner edge issue by discussing the behavior of six differently 
defined "inner edges" of slim accretion disks around a Kerr BH. We have found that the 
slim disk inner edges behave very differently than the corresponding Shakura & SUN- 
YAEV (1973) and NoviKOV & Thorne (1973) ones. The differences are quahtative 



74 



Relativistic, stationary slim accretion disks 





Figure 4.21: Profiles of tlie effective optical depth t^s for a = 0.01 and a* = in three 
different regimes of accretion rate. Vertical lines denote the locations of the BH horizon 
(dotted) and ISCO (dashed). Three types of behavior of can be seen: sharp drop at 
ISCO at the lowest accretion rates, maximum near ISCO for moderate accretion rates, 
and monotonic at all radii for the highest accretion rates. 
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and become important for the same range of luminosities independently of the BH spin. 
Even at moderate luminosities, L/L-Edd ^ 0.6, there is no unique inner edge. Differently 
defined edges are located at different places. For nearly Eddington luminosities, the 
differences are huge and the notion of the inner edge loses all practical significance. We 
summarize the properties and locations of the six inner edges in Table 4.4.7. 



4.5 Spectra 

The rate of viscous heating at given radius depends (in the a approximation) on the 
disk pressure and the rate of shear. These parameters change with radius what implies 
that the heating (and the corresponding cooling rate) significantly depends on radius 
too (profiles of the emitted flux of energy are presented e.g., in Fig. 4.2). In the stan- 
dard framework (assuming e.g., large optical thickness), the gas is assumed to radiate 
locally as a blackbody with the characteristic effective temperature corresponding to the 
amount of energy released in the form of the radiative flux. Since this quantity changes 
with radius, the outcoming spectrum of an accretion disk is a superposition of the black 
body spectra with a wide range of effective temperatures (multicolor blackbody). 

Let us consider the spectrum of a radiatively efficient disk. The outcoming flux of 
energy is given by the corresponding formulae of Shakura & SUNYAEV (1973) and 
NOVIKOV & Thorne (1973). For the non-relativistic case it is given by (Eq. 3.18), 




(4.46) 



n4 _ 

It is assumed that at given radius disk radiates the blackbody radiation corresponding 
to the effective temperature Tes- 

2h 

^^(^) = ^iW^K^- (4.47) 

The observed spectrum may be calculated by integrating the specific intensities (1^ = 
B„) over the whole surface of the disk*^: 

f cos ?' /"^o"' 

Su{r) = hdn = -— B,27rrdr, (4.48) 



disk 



where i is the inclination angle and D is the distance to the observer. 

Useful quantitative results may be obtained if we neglect the variations of Tcs near 
to the inner edge which have little impact on the shape of the total spectrum. Let us 
then approximate the radial profile of the effective temperature with 



> -p 
r ^ 



%S = Teff,in ( — ) , (4.49) 



^This formula neglects the gravitational g-factor. The precise way of calculating spectrum is given 
in Section 4.5.2 



^pot ^son '"var '"str '"rad '"rcf 

L/LEdd < 0.6 Tin ^ Tpot ~ r^on ~ '"var ~ '^str ~ ^^rad ~ '"rcf ~ '^ISCO 



L/L^dd > 0.6 



for a < 0.1 moves 
inward with 
increasing rh down to 

for a > 0.1 and 
sufficiently high rh 
disk enters the Bondi 
regime — undefined 



departs from 
ISCO; 

for a < 0.1: 

' son ~ ' mb) 

for a > 0.2: 

^son > nsCO 



undefined 



moves inward 
with increasing 
L down to BH 
horizon. 



moves inward 
with increasing 
L down to BH 
horizon. 



for 0.6 < 
L/LEdd < 1.0: 

J^ref < nsCO 

for 

L/LEdd > 1.0: 
undefined 



Table 4.1: Inner edges of a slim accretion disk: fpot - potential spout, Vson - sonic radius, r^ar - variability edge, Tgtr - stress edge. 
Trad - radiation edge, r^^i - reflection edge. 
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where Teg, in is the effective temperature at the inner edge (rjn = tisco for thin disks). 
According to Eq. 4.46 the exponent p for radiatively efficient dislcs is p = 3/4, while 
for advection-dominated flows (e.g., slim disks in the limit of high accretion rates) it is 
equal p = 1/2. Under this approximation, the spectral flux is given by. 



eff.in 



COS2 47r/irF„ //tsTeff.inV^^ 3 r°"'x(2/P)-l 



m 



]j2 p y hu 



u-^ / — dx, (4.50) 



with 



Xin = (4.51) 
.„.. = -4J-(^)'. (4.52) 

Kb ^cS;m \ Tin J 

The observed spectrum consists therefore of a power-law component 

oc z/3-(2/p) (4.53) 

limited by the low- and high-energy cutoffs: 

Mkiii f i^^y « ^ « Mkiii. (4.54) 

h \routJ h 

For a standard thin disk with p = 3/4, the spectrum is given by oc z/^^'^, while for 
advection dominated accretion (p = 1/2) it is oc . In Figure 4.22 we show an 
exemplary spectrum for p = 3/4 case. Its middle part (z/ = 10^^ ^ 10^^ Hz) follows the 
oc u^/'^ profile while the low and high energy tails follow the corresponding limits of 
the Planck function (Eq. 4.47) for the effective temperatures at Tout and r^n, respectively. 
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log V 



Figure 4.22: Typical spectrum of a thin accretion disk (p = 3/4) in a soft X-ray 
transient. Figure after Kato et AL. (2008). 
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4.5.1 Spectral hardening 

Real accretion disks do not radiate as a perfect blackbody. The emission is affected by 
such factors as the absorption opacity, density structure of the plasma, optical depth, 
rate of electron scattering, etc. Out of these effects most important in case of disks 
in X-ray binaries is the Compton scattering on hot electrons in the disk atmosphere. 
It distorts the spectrum towards temperatures that are higher than the effective tem- 
perature at given radius. However, when integrated over the disk, the total emergent 
disk spectrum can still be approximated to a good accuracy by a multicolor disc black- 
body, but with different characteristic temperatures (Shimura & Takahara, 1995; 
Merloni et al., 2000) 

T = /cTeff, (4.55) 

where /c is called the color temperature correction, or spectral hardening factor, since 
it is always greater then unity. In general, is a function of disk luminosity and its 
particular value depends on the disk structure close to its photosphere. 

To calculate the proper, hardened spectrum of an accretion disk one has to perform 
full radiative transfer calculations. DAVIS ET AL. (2005) calculated disk atmosphere 
models (BHSPEC) spanned on a grid of effective temperatures, effective gravity and 
surface densities relevant to disks in X-ray binaries. In the spectral model described in 
the following section, we either assume some particular value of the hardening factor 
/c or take the local spectra directly from the BHSPEC model. 

In Fig. 4.23 we present emission profiles of slim disk solutions for a* = 0.65 and three 
accretion rates: fn = 0.03, 0.3 and 0.6. For the highest accretion rate (corresponding to 
luminosity ~ l.OLEdd) the inward shift of emission caused by advection is visible. The 
dashed lines present profiles of the relativistic standard thin disk solutions (NOVIKOV & 
Thorne, 1973). The following figures (Figs. 4.24 and 4.25) present spectra calculated 
basing on these emission profiles. The spectra were obtained using ray-tracing routines 
developed by BURSA (2005) assuming either constant value of the hardening factor 
(/c = 1.7, Fig. 4.24) or using BHSPEC atmospheric models (Fig. 4.25). The slim 
disk solutions give in general slightly harder spectra, as they reach higher maximal 
effective temperatures in comparison to NOVIKOV & Thorne (1973) model (Fig. 4.23). 
Detailed discussion of the slim disk spectra will be given in BURSA & Sy\DOWSKl 
(2011). 

4.5.2 slimbb XSPEC fitting routine 

Using the numerical solutions of the slim accretion disk equations and an advanced ray- 
tracing technique (BURSA, 2005), we have calculated tables of spectra and developed 
a new model called slimbb for use with the X-ray spectral fitting package XSPEC 
(Arnaud, 1996). This model provides spectra of slim disks on a three-dimensional 
grid of the black-hole spin parameter, disk luminosity and inclination angle. They were 
calculated assuming either constant hardening factor fc, or using the BHSPEC disk 
atmosphere model (DAVIS ET AL., 2005; DAVIS & HUBENY, 2006). In the following 
subsections we describe most important properties of the model. 
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radius [GM/c^] 

Figure 4.23: Examples of local flux profiles emitted from the surface of slim disk (solid) 
and NOVIKOV & Thorne (1973) model (dashed lines) for BH spin a* = 0.65. Fluxes 
are plotted for three different mass accretion rates m = 0.03, 0.3 and 1.0. 



Disk model 

Slim disk solutions (Section 4.3) provide radial profiles of the following quantities: total 
amount of emitted flux F, radial velocity V , specific angular momentum surface 
density S and disk height cos Qh- They are used to construct and tabularize spectra of 
slim accretion disks at a grid spanned on accretion rate, BH spin and inclination angle. 

disk spectrum as observed by a distant observer 

The specific flux density of the disk radiation as observed by a distant observer at 
energy i?obs = hv is given by 

Fohs{v) = J Icm{i^/g,r,fic)9^dQ , (4.56) 

where g is the photon redshift (Eq. 4.60) and dQ is the element of the solid angle 
subtended by the image of the disk on the observer's detector plane, /cm is the specific 
intensity of the radiation coming out from a point on the disk in the right direction to 
reach the observer as measured in the fluid comoving frame. Generally, it is a function 
of energy, place of emission and emission angle, slimbb model assumes for /cm either a 
diluted blackbody or takes the emerging radiation spectrum from BHSPEC. 

In the case of a diluted blackbody radiation, the specific intensity is given by Planck's 
law with color-modified temperature 

= , , T(/ie) . (4.57) 

The effective temperature T^s can be obtained from the total energy flux F: 

4 

^cff = I I . (4.58) 



out 

cr 
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A/ = 10A/o, a = 0.65 
i = 0.1, inc = fi6° 
hf = 1.7 




energy [keV] 



Figure 4.24: Calculated theoretical continuum spectra corresponding to the emission 
profiles presented in Fig. 4.23. The plot shows spectra of slim disk (solid) and NOVIKOV 
& Thorne (1973) (dashed lines) models calculated assuming constant value of the 
hardening factor (/ = 1.7) for an observer at inclination 66°, and distance lOkpc. 



The color hardening factor fc accounts for the effect of free electrons in the disk 
atmosphere on photons and the fact that the locally observed color temperature of the 
emitted radiation is generally higher than the effective temperature of the disk related 
to the locally dissipated energy (see Section 4.5.1). T(yUe) describes the dependence on 
the photon emission angle (the angle between the emitted photon wavevector and the 
disk surface normal), which is T = l for isotropic emission or 

T = ^(2 + 3/ie) (4.59) 
if the emission is assumed to be limb-darkened. /Xg is the cosine of the emission angle. 



g-factor and disk emision angle 

Most of the flux emitted by disk comes from the region of strong gravity and high 
velocities, and is subject to intense Doppler boosting and gravitational redshift. 

The combined effect of gravitational and Doppler energy shift is described by g- 
factor, which is the ratio between energy at which the photon is observed and energy 
at which it has been emitted. It is calculated using the general formula. 



where p'^ is photon four-momentum and we obtain u'^ by a transformation from the 
co-rotating frame, = e|*^^ u^'^\ The fluid velocity measured by a co-rotating observer 
has only the radial component V and we have 



= (Vl + y^ -V, 0, 0) 



(4.61) 
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Figure 4.25: Same as Fig. 4.24 but calculated using the BHSPEC atmospheric model 
instead of assuming constant value of the hardening factor. 



Here, V, although negative, is taken with the minus sign as the radial vector e^^^ of the 
orthonormal frame points in the inward direction. 

The emission angle is the angle between the direction of the normal to the disk 
surface and the direction of escaping photon momentum. The cosine of the emission 
angle is given by the general formula 

/.e = , (4.62) 

where n'^ is the normal vector to the disk surface. The disk normal can be easily 
expressed in the co-rotating frame, where it is 

^(a) = g(a)^M = (0, sin 7, COS 7, 0) . (4.63) 

The angle 7 between the disk surface normal and the direction of the poloidal basis 
vector e^^-* is 

7 = ^-(^ + «), (4.64) 

where 'd is the poloidal coordinate angle of the place of emission and a is the disk height 
derivative with respect to the cylindrical radius r, 

a = tan-i . (4.65) 

In Boyer-Lindquist frame, the disk normal vector, normalized to unity, is then 

= [0, cos i^ + a), ^— sin (^9 + a), ) . (4.66) 



An example of how slim disk changes the profiles of (7-factor and emission angle is 
given in Fig. 4.26. For a chosen observer's viewing angle of 66° and BH spin a* = 0.6 
we compare these two quantities for two different mass accretion rates. 




impact a [GM/c.^] impact a [GM/c^] impact a [GM/c^] 

Figure 4.26: Profiles of emission angle (top row) and (7-factor {bottom row) for two mass accretion rates m = 0.3 and 0.6. Other 
parameters like spin a* = 0.6, alpha viscosity a = 0.1 and inclination i = 66° are fixed. The two leftmost columns represent slim 
disk, while the right column show for comparison corresponding thin disk picture (NT solution), which looks the same for any 
accretion rate. For slim disk the mass accretion rates correspond to luminosities O.SLEdd and l.OZ/Edd! respectively. Plots show 
the image plane view of the disk. 
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For luminosities less than O.SLEdd (mass accretion rates M < 0.2 for a* = 0.6) the 
two models differ only slightly. For higher luminosities, presented in Fig. 4.26, the 
increasing difference between these models comes from the fact that the radiation inner 
edge of slim disks shifts to smaller radii, that the disk thickness increases, and due to 
more and more significant departure of angular momentum from the Keplerian profile 
(see Section 4.3). 

Self-obscuration 

We may notice an interesting consequence, unknown to thin disks, namely that at this 
inclination the inner regions of slim disks may be obscured by its outer parts. This 
fact is noticeable in the plot of the emission angle, where we see regions from which 
the radiation comes almost parallel to the disk surface. Full obscuration of the central 
BH takes place only for super-Eddington luminosities. This effect has also obvious and 
important implications for the temperature-luminosity relations. 

The condition for self-obscuration is roughly that h/r somewhere in the disk is larger 
than the cosine of inclination. If for higher luminosities the radiation pressure becomes 
important, the h/r ratio typically has a maximum near radius where radiation pressure 
dominates over gas pressure. The most luminous part of the disk, however, lies still 
inside this radius and thus if the condition H/R> cost is met, the strongly radiating 
part of the disk disk becomes partially or even fully obscured by an outer part. This, 
of course, has important observational consequences. 

Let us consider our model example of a source with disk inclination 66°. The critical 
h/r ratio is 0.4 and the situation is illustrated in Figure 4.27. The plot shows the profiles 
of the disk photosphere for several luminosities and compares them with trajectories of 
photons sent from infinity with zero azimuthal angular momentum. For zero spin case 
(top panel), we see that already at luminosity O.SLgdd photons from the part of the 
disk near the BH need to escape almost parallel to disk surface to reach the observer. 
For even larger luminosities this part of the disk starts to be completely obscured and 
at 2 -^Edd the large part of the inner disk is eclipsed including the BH itself. For high 
spin case (bottom panel), the situation is quite similar, only the luminosity has to be 
higher. 

Model parameters 

The SLIMBB model comes in a form of a FITS table with precalculated spectra and a 
small routine which reads the data and produces the final spectrum in terms of linear 
interpolation in the parameter space. For five different values of a viscosity the table 
contains extensive three dimensional grids of spectra spanned on the spin, luminosity 
and inclination. These parameters vary within the limits given in Table 4.2 with steps 
chosen to ensure that the spectra differ roughly equally between any two adjacent steps. 

The model has all together nine parameters with the following meaning (see Ta- 
ble 4.2 for their limit values): 

M - mass of the central BH in units of solar masses, 

a - specific angular momentum of the central BH in units of GM/ c. 
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Figure 4.27: Illustration of disk self-obscuration in the z-x plane. The figure shows 
vertical cuts through the accretion disk, location of disk photosphere for different lumi- 
nosities and trajectories of photons with zero azimuthal angular momentum at infinity 
reaching an observer at inclination 66deg. Top panel corresponds to spin a* = 0.0, 
bottom panel to spin a* = 0.9. 

L - disk luminosity in units of Eddington luminosity L-^dd, 

inc - inclination of the observer with respect to the symmetry axis in degrees, 

D - distance of the source in kiloparsecs, 

f_hard - hardening factor (if positive) or a switch for BHSPEC (if negative), 

alpha - the value of a-viscosity parameter, 

Iflag - switch for isotropic (0) or limb-darkened (1) emission, 

vflag - switch for equatorial plane (0) or photospheric (1) ray-tracing. 
Most up-to-date version of SLIMBB and its detailed description may be found at 
http : / /astro . cas . cz/bursa/pro j ects/slimbb. 
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Table 4.2: 


slimbb model parameters. 




Parameter 


Units 


Default value 


Min value 


Max value 


M 


Me 


10. 


0. 


100. 


a 


GM/c 


0. 


0. 


0.999 


L 


-^Edd 


0.1 


0.0003 


1.0 


inc 


deg 


0.0 


0.0 


85.0 


D 


kpc 


10. 


0. 


1.0e4 


f hard 




1. 


-10. 


+10. 


alpha 




0.1 


0.005 


0.1 


Iflag 




0. 


-1. 


+1. 


vflag 




0. 


-1. 


+1. 
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Chapter 5 



Relativistic, non-stationary slim 
accretion disks 



X-ray binaries and AGN exhibit different types of variability taking place at different 
timescales. Some of them may result from intrinsic instabilities of accretion disks. 
Other may result e.g., from the orbital motion. Instabilities caused by the viscous 
and non-adiabatic processes in disks may be classified into three types (Kato ET AL., 
2008): thermal, secular instabilities and overstabilities of waves and oscillations. The 
first one takes place on the thermal timescale what makes it most profound. Its origins 
will be discussed in the following paragraph, while in the next sections we will present 
numerical model of non-stationary slim accretion disks which is able to reproduce the 
limit cycles resulting from the thermal instability. 



5.1 Thermal instability 

Thermal instability takes place when a temperature perturbation over the equilibrium 
is amplified by a thermal process. Therefore, the characteristic timescale is the thermal 
timescale (Kato et al., 2008), 

^2 ^2 

n. = (5.1) 

over which the disk temperature changes. The following relations are satisfied, 

h 

T"hyd = — < Tth < r^is = — , (5.2) 

where Thyd is the hydrodynamical timescale, over which the vertical disk structure varies, 
while Tvis is the viscous timescale determining the radial disk structure changes. 

When studying the secular and thermal instabilities one has to filter out the acoustic 
oscillations occurring on timescales shorter than both thermal and viscous ones. It 
can be done by neglecting the time derivatives in the momentum equations. On the 
contrary, as the instability involves temperature perturbations, the energy equation 
should be taken into account in a time-dependent form. The dynamical equilibrium 
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may be adopted as it is achieved over the dynamical timescale, much shorter than the 
thermal one. Finally, since we consider the force balance between the gravitational and 
centrifugal forces, we may assume the gas motion resulting from a temperature change 
occurs only in the vertical direction (S = const). 

Under these assumptions we may perform a simple thermal stability study. Full 
mathematical derivation of the stability criterion for both thermal and secular instabili- 
ties has been given e.g., in Shakura & SUNYAEV (1976) and PiRAN (1978). Recently, 
ClESlELSKi ET AL. (2011) have generalized the perturbative analysis of the stability 
criterion to account for modified (delayed, advanced or with different stress-to-pressure 
response factors) viscosity prescriptions. 

Let us write the vertical equilibrium equation in the following form (compare Eq. 3.13), 

p = ^Q^m. (5.3) 

In the innermost region of an accretion disk we have p ^ Prad and n ^ n^s = const. The 
viscous heating rate is given by, 

Qvis ^ _^Qtr^h_ (5.4) 

Under the ap assumption (Section 2.2.1) and taking Eq. 5.3 into account, we have, 

g"" (xphoc h\ (5.5) 
The radiative cooling rate per unit surface is, 

grad _ ^ (5 g) 

where we used again Eq. 5.3. If the temperature exceeds slightly the equilibrium (g"*^ = 
Q^^'^) value, the disk expands in the vertical direction as /i oc p oc (Eq. 5.3). As a 
result, both Q^^^ and Q'^^'^ increase, but the rate is larger for the heating component. 
This means that the specific entropy of the gas (and temperature) increases, which 
leads to a further increase of H. In this way, the positive feedback occurs and the disk 
becomes unstable. 

In the gas dominated case, the heating rate is still given by, 

g^'^ oc h\ (5.7) 

However, the radiative cooling is now expressed as, 

Q'""^ oc /^^ (5.8) 

since pj-ad oc T'', p = pg^s oc Tp and hence prad oc {p/pY oc {phY oc . For such 
dependence of the heating and cooling on if, the reasoning given above will not be 
valid. The disk in the gas pressure-dominated case is thermally stable. 

The general criterion for the thermal stability may be put in the form (Pringle, 
1976), 

|^(gvis_grad 



dh 



> 0. (5.9) 
s 
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Using Eq. 5.3 we can relate the smaU-amphtude perturbation of pressure pi superposed 
on the equihbrium value Pq as, 

— = — , 5.10 

Po ho 

where similar perturbations of the disk thickness have been introduced. Using the 
general equation of state (Eq. 3.21), we may also write, 

Ti _ 1 + 13 h 

To "4- 3/3 /.o' ^^-''^ 

where /3 = Pgas/p is the gas pressure to total pressure ratio. The heating rate, Q"'^^, 
follows Q^"* oc ph and therefore, 

^ = 2— (5.12) 

Qt ho 

Similarly, the perturbation of Q^^'^ follows Q^^'^ oc and thus, 

Q-'^ _ 4(1 + /3) h 
Qr~ 4-3/3 /.o- 

Putting the equations given above into Eq. 5.9 we obtain the stability criterion, 

2-5/3<0. (5.14) 

The regions of accretion disks with /3 < 0.4 are thermally unstable. This is the case if 
only the accretion rate exceeds some critical value, e.g. m 0.06 for a non-rotating 
BH. 

Fig. 5.1 summarizes results of a more detailed, based on the perturbative analysis, 
study of thermal and secular types of instabilities in accretion disks. The growth rates 
of perturbations (solutions of the dispersion relation) for some arbitrary values of /3 are 
plotted against the perturbation wave number k = 27t/X. Positive values correspond to 
the unstable case, while negative prove stability. It is clear that only for /3 < 0.4 the 
solutions lie above the x-axis. These curves have two branches: the lower (approaching 
zero for A — t- oo) corresponds to the secular instability, while the upper to the thermal 
one. For /3 = 0.4 the disk is marginally stable in the long wavelength limit and stable 
for finite wavelengths. Higher values of /3 give solutions with the negative sign reflecting 
their stability. 

The stability criterion (Eq. 5.14) is very strong. It predicts that only disks with very 
low accretion rates are stable. One could therefore expect strong variability of this type 
in most accretion disks. However, it is not the case since most of observed microquasars 
exhibit in their high/soft states stable emission for periods of time much longer than the 
corresponding thermal timescales. Moreover, most recent MHD simulations (HiROSE 
ET AL., 2009a) show that thermal instability may not occur in real accretion disks. 
This inconsistency has not yet been resolved. Recently, ClESIELSKI ET AL. (2011) have 
addressed this issue, showing that time delays in heating or reduced stress-to-pressure 
response may quench the thermal instability. 

Considering the unique limit-cycle behavior of GRS 1915+105, which contains near- 
extreme Kerr BH and exhibits limit-cycle-like luminosity variations, we suggest that the 
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Figure 5.1: Growth rate as a function of the radial wave-length of perturbations for 
a = 0.1, 7 = 5/3 and different values of /3. The branch closer to the horizontal axis 
denotes the secular mode, while the other represents the thermal one. Dotted curves 
denote the complex conjugate solutions of the dispersion relation. 

classical theory may be indeed correct in predicting the radiation pressure instability, 
but under some, so far unknown, conditions this behavior is suppressed, as observed 
in most X-ray binaries. Since the proper description of viscosity is still developing, the 
limited classical theory is the only reasonable choice for our study before a widely ac- 
cepted new theory is presented. In the next sections we investigate the time-dependent 
behavior of relativistic slim accretion disks that exhibit limit cycles resulting from the 
classical thermal instability. 



5.2 Equations 

In this Chapter, we consider axisymmetrical relativistic accretion flows onto Kerr BHs. 
We use the Boyer-Lindquist spherical coordinates (t, r, 6, (p) to describe the spacetime 
around a BH. Putting all of the complicated derivations into Appendix B, we give here 
the basic equations governing the dynamical behavior of the flow. Assuming dt = 0, 
they all reduce, upto the coefficients related to the vertical integration, to the equations 
given in Section 4.1. 
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(5.20) 



where A = — 2Mr + a^, A = + r^a^ + 2Mra^ and 7 is the Lorentz factor. 
Equations 5.15, 5.16, 5.17, and 5.20, defining the time derivatives of the surface density, 
radial velocity, angular momentum and central temperature, respectively, represent 
corresponding conservation laws. Equations 5.18 and 5.19 determine the evolution 
of the disk thickness and the vertical acceleration. We adopt the diffusive form of 
viscosity (Eq. 2.20) in equations 5.17 and 5.20. In these equations, surface density S, 
radial velocity V , specific angular momentum £, half thickness /i, vertical velocity U , 
and temperature T, are six essential quantities describing the structure of disk. 
The other symbols denote the following auxiliary quantities. 
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where equation (5.26) gives a bridging formula that is valid for both optically thick and 
thin regimes (HuBENY, 1990), a is the viscosity parameter, prad is the radiation pres- 
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sure, and r^f and Tp are the Rosseland and Planck mean optical depths (SZUSZKIEWICZ 
& Miller, 1998). 

The properties of an accretion flow depend on three crucial parameters: 



Dimensionless BH spin a*: 



a. = ^, < a, < 1; (5.28) 



• Diffusive viscosity parameter a (Eq. 5.21): 

< a < 1; (5.29) 

• Dimensionless mass supply rate m (compare Eq. 4.35): 

m = — : . 5.30 

MEdd 

Where M is the BH mass and M(rout) is the accretion rate at the outer boundary 
{r = rout). 



5.3 Numerical methods 

The given set of equations may be put in the following, general form, 

dui{r, t) 



dt 



Huk{r,t)) (5.31) 



where Ui{r,t) stands for all six functions involved (S, V, C, H, U and T) and L is 
a partial differential operator. Such a set of equations may be solved numerically by 
discretizing in time and radius. For the purpose of solving the non-stationary equations 
of slim disks we applied the standard Chebyshev pseudospectral method for the spatial, 
and combined Runge-Kutta and third-order backward differentiation explicit schemes 
for the time discretization, respectively. 

5.3.1 Spatial discretization 

A physical quantity u{r) is approximated by the following series with a finite number 
of terms: 

U{rk) = U{g{fk)) = ^n=0^nTn{fk) = ^n=0^n COS (^^^^ , (5-32) 

where Tn{fk) is the n^'^-order Chebyshev polynomial; = cos^kir/N) [k = 0, 1, ■■.,N) 
represents the Chebyshev- Gauss-Lobatto collocation points, with being the number 
of collocation points; = g{fk) is the mapping from the Chebyshev-Gauss-Lobatto 
collocation points — 1 < "Tfe < 1 to the physical collocation points rmin "^"^k ^ ''"max and 
is a strictly increasing function that satisfies both g{—l) = and 5^(1) = r^a.^; and 
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the Un are the spectral coefficients and can be calculated from the physical values u{rk) 
with a fast discrete cosine transform (FDCT; PRESS (2002)). 

In a similar way one may approximate radial derivatives, 
du{rk) 1 du{g{fk)) 1 ^/ 



dr dg/dr dr dg/dr 



Silo<^n(rfc), (5.33) 



where the spectral coefficients u'^ are evaluated basing on u„ using the following recur- 
sive relations. 



u'^_^ = 2Num, (5.34) 

= M^+2 + 2(n + l)Mn+l, 



with Co = 2 and c„ = 2 for n = 1,2, N. Calculation of these coefficients is followed by 
an inverse FDCT of providing du{g{fk)) / df and, finally, discrete spatial derivatives 
du{rk)/dr. 



5.3.2 Time discretization 

The integration in time is performed using two algorithms. First two time-steps are 
integrated using the third order total variation diminishing (TVD) Runge-Kutta scheme 
(Shu & OSHER, 1988), while the further time-evolution using the third order backward- 
differentiation explicit scheme (Peyret, 2002) which provides good enough accuracy 
and is less CPU-consuming. 

The third order TVD Runge-Kutta scheme may be expressed as, 

+ AtL(M"), 

^w« + lw(i) + lAtZ(w(i)), (5.35) 
in" + ^n(^) + ^AtL(n(^)), 

where At is the time-step; and m""*"^ are the values of the physical quantity u at 
the n-ih. and (n + l)-th time-levels, respectively; and m^^^^ and u^"^^ are two temporary 
variables. 

The third order backward-differentiation explicit scheme is defined as follows, 

a,u-^'~^ = bMu--'), (5.36) 

j=0 j=0 
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(5.38) 



with t", and being the values of time of the n-th, (n — l)-th, and (n — 2)-th 
time-steps, respectively. This scheme consumes about three times less computational 
time than the Runge-Kutta scheme but is not able to start the time-integration by 
itself. Therefore, we combine these approaches and start integration with the latter, 
followed by the faster third order backward-differentiation explicit scheme. In such a 
way we are able to obtain a sufficient high order accuracy with minimal CPU-time 
consumption. 



5.3.3 Specific techniques and assumptions 

Besides discretizing the equations describing non-stationary evolution of an accretion 
disk (Eqs. 5.15 - 5.20) one has also to treat properly several numerical issues which may 
affect the integration. In order to obtain a physical, transonic and numerically stable 
solution in a finite domain, it is necessary to impose appropriate boundary conditions 
and to apply some filtering techniques to overcome the inevitable spurious nonlinear 
numerical instabilities in the code. For the purpose of this project we applied the 
following techniques: 

(i) domain decomposition into 7 subdomains so that each sub-domain contains at 
most one single region of rapid variation, 

(ii) spectral filtering by using an exponential filer in space to filter out the high- 
frequency modes in each time step, 

(iii) small numerical viscosity which makes the lowest density regions stable. Details 
of these algorithms are given in Ll ET AL. (2007). 
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In this section we briefly discuss the hmit-cycle evolution of an accretion disk for a 
non-rotating BH and mass accretion rate high enough to trigger the instability (we 
choose rh = 0.06) assuming a = 0.1. Such a disk is radiation pressure dominated in its 
inner parts and is expected to exhibit the thermal instability. 

The top-most panel of Figure 5.2 presents the disk thickness just before the start 
of the cycle (t = Os). As the accretion rate is relatively low, the disk is geometrically 
thin {h/r ^ 1). The disk structure does not correspond precisely to the stationary 
solution of a slim disk with the same model parameters, as the fluid configuration is 
not stationary and the disk is in the process of reformation which is terminated by 
the onset of the instability. As it sets in {t = 2s, the second panel of Figure 5.2 and 
corresponding lines in Figure 5.3) the temperature rises rapidly, the disk expands in the 
vertical direction and a shock appears in the surface density and accretion rate profiles. 
These shocks move outward with time, forming an expansion wave, heating the inner 
material, pushing it into the BH, and perturbing the outer gas. It is worth mentioning 
that the accretion rate at the shocked region (the border between the expanding medium 
and the quasi-stationary disk) is negative which means radial outflows, causing the 
outer, incoming gas to further pile up. 

At t = 12s (e.g., the third panel of Fig. 5.2) the disk thickness and the negative 
spike of the mass accretion rate in the expanding region reach their maximum values. 
The local rh inside the wavefront, however, exceeds 3 which is far above the initial value 
m = 0.06 and reflects the rapid inflow of gas. 

Once the wavefront has moved beyond the unstable region (r < 120rg = 240M for 
this case) the disk deflates as the instability is no longer supported by radiation pressure 
dominated parts of the disk. However, the expanded region still moves outward. The 
temperature and surface density in the inner parts rapidly drop down. The inner region 
becomes geometrically thin, with the temperature and surface density lower than the 
corresponding values at t = Os 

The outburst is essentially over after the thermal timescale (e.g., t = 32s). It is 
followed by a much slower process (taking place on the viscous timescale) of refilling 
and reheating of the inner part of the disk. Finally, after ~ 976s, the disk returns to 
the same state as at t = Os and the thermal instability triggers again leading to a new 
cycle. 

The bottom-right panel of Fig. 5.3 presents the bolometric luminosity of the disk, 
obtained by integrating the radiated flux per unit area over the disk at successive times 
for three complete cycles. The light curve exhibits a burst with a duration of about 20 
seconds followed by a quiescent phase corresponding to the disk reformation process. 
During the outburst the disk may reach super-Eddington luminosities (1.5LEdd for this 
model). 
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Figure 5.2: Evolution of the disk thickness during one full cycle for m = 0.06, a = 0.1 
and a* = 0; = 2M. Figure after Ll ET AL. (2007). 



5.4 Limit cycles 



97 





Figure 5.3: Evolution of the temperature (top-left), surface density (top-right) and 
radial velocity (bottom-left panel) during one cycle. The bottom-right panel presents 
the outcoming light curve covering 3 cycle durations. A non-rotating BH, m = 0.06 
and a = 0.1 were assumed. Figure after Ll ET AL. (2007). Vg = 2M. 
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Time 



Figure 5.4: Sketch for explaining the definitions of characteristic quantities. 

5.5 Results 

5.5.1 Parameter space and characteristic quantities 

In order to understand the effects of BH spin on the hmit-cycle behavior, it is necessary 
to investigate the parameter space (a*,Q;,m) with a series of models. In this section, 
we describe and discuss results of twelve runs assuming different sets of the input 
parameters span on the grid: a* = (0,0.5,0.95), a = (0.07,0.1), and m = (0.06,0.1) 
(cf. Table 5.1). All of the cases are thermally unstable and undergo recursive limit-cycle 
evolution. 

To describe the limit-cycle behavior, we define four characteristic quantities, which 
are the Outburst Duration (the full width at half maximum of light-curve, hereafter 
OD), Cycle Duration (time interval between two outbursts, hereafter CD), ratio of OD 
to CD, and Maximal Luminosity (the peak value of the light curve during outburst, 
hereafter ML; calculated intrinsically, with no ray-tracing), respectively. In Figure 5.4, 
we visualize these definitions on a sketched light-curve. 

The computations for each one of those twelve cases were continued until the CD 
converged to a constant value. The constancy of CDs implies that the computations 
have reached the state uniquely determined by the crucial parameters (a*,Q!,m), with- 
out any impact of the initial conditions. In practice, the required level of convergence is 
achieved after three or five cycles. We performed computations for each case for more 
than seven cycles (some cases reached ten cycles) and we based our subsequent analysis 
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Table 5.1: Limit-cycle characteristic quantities. See Section 5.5.1 for the definitions of OD, CD, and ML. The values showed 
the columns, OD, CD, OD/CD and ML are collected as "mean value" ± "standard deviation" (relative deviation). 
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on the last four cycles only. 

In Table 5.1 we list the mean values, standard and relative deviations of each char- 
acteristic quantity for our twelve models. These values were calculated basing on the 
retained cycles only. Since the relative deviations (the percentages in the brackets) are 
all very low, we do not show the error bars for the data points in the subsequent figures. 



5.5.2 Impact of BH spin 

In order to reveal the effects of BH spin on limit-cycle behavior, we collect our twelve 
models into four groups to reflect the different effects of viscosity and mass supply (i.e., 
mass accretion rate at the outer boundary). The models with (a = 0.07, m = 0.06) 
and different spins are represented by empty stars, (a = 0.1,m = 0.06) by filled stars, 
{a = 0.07, m = 0.1) by empty circles, and {a = 0.1, m = 0.1) by filled circles in the 
Figures 5.5 and 5.6 (both are based on the data from Table 5.1). 

In the four panels of Figure 5.5 we plot OD, CD, OD/CD and ML versus a*. It is 
clear that there is no distinct and consistent dependence on neither for OD nor CD 
(panels (a) and (b)). The ratio of OD to CD (panel (c)) is almost independent of for 
all sets of input parameters. On the contrary, ML exhibits a perfect correlation with 
a* (panel (d)), though this relation is still slightly distorted by the impact of a and rh. 

These facts may be easily understood in the framework of the general relativity. The 
effects of BH spin are restricted to the innermost parts of accretion disks. Therefore, 
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Figure 5.6: Plot of maximal luminosity versus ratio of OD to CD. The numbers near 
each data point are the values of a*. The meanings of markers are the same as those 
in Figure 5.5. The dashed lines are the OD/CD-ML fitted lines for different a* cases. 
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the quantity ML, which corresponds mostly to the radiation emerging from the inner 
part of an accretion disk, must display a strong dependence on a*. The other quantities 
(OD, CD and OD/CD) cannot have such strong dependence on a* as they depend on 
the disk structure and its evolution at wide range of radii (from the inner part up to 
more than 200M), where the effects of viscosity and mass-supply dominates over the 
effect of BH spin. 

5.5.3 Possibility of estimating the black hole spin with limit- 



As has been discussed in the previous section, the effects of BH spin on the quantities 
OD, CD and OD/CD are easily disturbed and even obscured by the complex effects of 
viscosity and mass-supply. Thus, they cannot serve as proper probers of BH spin. Even 
ML, which has monotonic and positive correlation with a*, cannot be used directly to 
estimate the BH spin, because the ML dependence on BH spin is significantly affected 
by other factors (see panel (d) of Figure 5.5). 

However, it is reasonable to base both on ML and the OD to CD ratio. In Figure 
5.6, we show the OD/CD-ML diagram for all the models. The OD/CD ratio raises with 
increasing values a and m and is hardly dependent on a* (see panel (c) of Figure 5.5). 
It is remarkable, that the combined effects of viscosity and mass-supply make ML scale 
almost linearly with OD/CD ratio. For a given BH spin there is a single straight line 
reflecting the ML dependence on OD/CD ratio. These lines do not intercept and have 
similar slopes. The OD/CD-ML diagram may be, therefore, used for rough BH spin 
estimation if only the limit cycle parameters are read from observational data. 

We find that the dependence of ML on OD/CD may be approximated by the fol- 
lowing formula 



where f] is the efficiency of accretion for a thin disk and is given by (Eq. 2.29), 



with risco being the radius of the marginally stable orbit. In Figure 5.6 we plot with 
dashed lines the fits obtained with these formulae for a* = 0, 0.5 and 0.95. Knowing 
the values of ML and OD/CD, one can easily obtain from Eq. 5.39 the radius of the 
marginally stable orbit and, subsequently, the rough estimate of BH spin. 

There are, however, at least few limitations for the application of this method. The 
first one is related to the fact that the OD/CD-ML relation is still subject to significant 
dispersion caused by the combined effects of a and m (see Fig. 5.6) as well as to 
ray-tracing effects we neglect, related to e.g., the inclination angle and limb darkening. 
Another one is the model-dependence of our method. We assume the limited classical 
theory to describe the accretion flow. We neglect mass outflows, magnetic fields and 
apply the a prescription of viscosity with a independent of radius. These factors, 
especially the viscosity treatment, may have significant impact on the ML vs OD/CD 
relation that we have obtained. Nayakshin et AL. (2000) showed that reproducing 



cycle 



ML/LEdd = 7.59 OD/CD + 0.71 + 7.18 77, 



(5.39) 




(5.40) 
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the light curves of GRS 1915+105 is possible only if much more complicated models are 
considered. Nevertheless, if the limit-cycle behavior is the actual reason for the observed 
variability of some microquasars, and if the extraction of ML and OD / CD parameters 
from the light curves is possible, then, under the limitations of the model presented 
in this work, the BH spin can be roughly estimated using the method presented here, 
before other accurate methods, e.g., fitting spectral energy distribution, are applied. 

5.5.4 Evolution of the temperature-luminosity relation for the 
limit-cycle in classical theory 

The maximal disk temperature-luminosity (T-L) diagrams have been used to compare 
the observations with the predictions of the classical theory of limit-cycles (e.g., GlER- 
LiNSKi & Done (2004); Kubota & Makishima (2004); Kubota et al. (2001)). 
In Figure 5.7, we show the T-L evolution for the last calculated cycle of each model. 
The panels present the evolution for a given set of a and m parameters and different 
BH spins {a^ = 0, 0.5 and 0.95 marked by solid, dashed and dotted lines, respec- 
tively). These curves reveal several common features: (1) All of the cycles begin at the 
closely-located triangle points and spend a quarter of CD evolving along the curve in 
the anti-clockwise direction before they reach the square points, and finally spend the 
other 3/4 of CD to return back to the starting location (triangles); (2) The evolution 
of all the cycles fits (with exception to the outbursts of the a* = 0.95 case) between 
two straight lines corresponding to Ljisk oc T^^^x relations with different color correction 
factors (/col = 1 and 7.2, see GlERLINSKI & Done (2004)); (3) The T-L curves of 
all the cycles are similar for most of CD except for the outburst state — the BH spin 
impact is revealed at the outburst of limit-cycle only; (4) During the mass restoring 
process (prior to the next outburst) the disk obeys L^isk oc T^l^ (see the dashed straight 
lines in all panels). 



5.6 Discussion 

The theory based on the standard viscosity prescription 

W = -aP, (5.41) 

where W and P are vertically integrated stress and total pressure, respectively, pre- 
dicts that radiation pressure-supported regions of radiatively-efficient accretion disks 
are thermally and viscously unstable. The unstable luminosity range is smaller for an 
(ad hoc) alternative viscosity prescription 

W = -a^/p^^^s, (5.42) 

where < /i < 1 is a constant. The value fi = corresponds to the standard pre- 
scription, and /i = 1 to the often considered "geometrical mean" prescription. The 
instability disappears for /i > 8/7. From the recent important paper of HiROSE ET AL. 
(2009b) it follows that the standard Shakura & Sunyaev (1973) viscosity prescrip- 
tion (/i = 0) fits much better results of their MHD simulations of accretion disks than 
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Figure 5.7: Evolution of the last cycle for the groups with different ot and m in the 
maximal disk temperature - luminosity (T-L) plane. All cycles begin at the triangle 
points, and then evolve along the curve in the anti-clockwise direction. The square 
points are used to denote the quarter of cycle durations. The solid straight lines in 
each panel represent L^isk oc Tj^^x ^ non-spinning BH (see Eq. (3) of GlERLINSKI 
& Done (2004)). They are calculated with the color temperature correction /coi = 1 
for the left lines and /coi = 7.2 for the right ones. The dashed straight lines denote the 
Ldisk oc Tj^'Jx relation. In each panel, the solid, dashed, and dotted curves are for the 
different BH spins with = 0, 0.5, 0.95, respectively. 
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the alternative prescriptions (5.42) with n > 0. However, HiROSE ET AL. (2009a) also 
found that the radiation pressure supported MHD disks were thermally stable in their 
simulations (it is not known whether they are viscously stable). This result contradicts 
the hydrodynamical stability analysis, and therefore needs an explanation — what are 
the MHD effects that stabilize the disk? HiROSE ET AL. (2009a) argue that although 
the viscosity prescription has the form (5.41), it should nevertheless be modified in a 
subtle way, because the "viscous heating" Q"*^ that occurs in the disk is delayed with 
respect to the MHD stress W = Wmhd, 

g--(t) = W(t-At)^. (5.43) 

ar 

Here At is the delay, found to be about 10 dynamical times, and Q is the angular 
rotation of matter. It should be obvious that this modification is irrelevant in the 
stationary case. 

In an exhaustive Newtonian analytic stability analysis (ClESIELSKI ET AL., 2011) 
our group confirmed (in general) the suggestion made by HiROSE ET AL. (2009a) that 
the delayed heating (5.43) stabilizes, in some parameter range, the radiation pressure 
thermal instability. However, we also found a variety of different parameter ranges with 
interesting, and complex, oscillatory behaviors that need to be further investigated. 

On the observational front, spins of BHs in several galactic BH binaries have been 
recently measured by fitting their spectral energy distribution to the relativistic geo- 
metrically thin accretion disk model (e.g., Shafee ET AL. (2006); McClintock ET 
AL. (2006)). From among these objects, the GRS 1915+105 is the most special one. It 
is believed to have a near-extreme Kerr BH with the spin very close the theoretical max- 
imal value a* = 1 (see, McClintock ET AL. (2006)). It is also the only object that 
exhibits the quasi-regular luminosity variations (Belloni ET AL., 1997; Nayakshin 
ET AL., 2000; Janiuk ET AL., 2002; Watarai & MiNESHIGE, 2003; Ohsuga, 2006; 
Kawata ET AL., 2006; Janiuk & Czerny , 2011), similar to the hmit-cycle predicted 
by the standard, Newtonian, hydrodynamical models. There are obvious questions to 
investigate here. May one explain the luminosity variations observed in GRS 1915+105 
as the classic "slim disk" limit-cycle behavior? What is the role of the nearly extreme 
spin of GRS 1915+105 in the context of the observed variability? Could the case of 
GRS 1915+105 help to find the answer about the proper viscosity prescription? 

Our research shows that the shape of light curves is weakly dependent on BH spin. 
The only parameter that feels the impact of the BH angular momentum is the maximal 
luminosity. Therefore, under the assumptions adopted in this paper (e.g., constant a 
in radius), limit cycles produce light curves only qualitatively similar to those observed 
in GRS1915+105, independently of BH spin. To obtain better agreement one has to 
construct more complicated models (e.g., Nayakshin et AL. (2000)), but the mech- 
anism of the instability may remain intact. We prove that the cycle duration weakly 
depends on the rotation of BH. Thus, the observed value for GRS1915+105 most prob- 
ably results from combined effects of viscosity and accretion rate and is not affected 
by the spin of the central BH. The very presence of limit cycle variability only in this 
object challenges our understanding of the mechanisms of viscosity. 
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Chapter 6 

Vertical structure of slim accretion 
disks 

Slim disks, just like Shakura & Sunyaev (1973) model, are based on the assumption 
that the dissipation mechanisms operating in accretion flows may be described by a 
viscous stress tensor whose leading component is proportional to the pressure. 

Slim disk solutions, again just like the Shakura & Sunyaev (1973) and NoviKOV 
& Thorne (1973) solutions, are obtained by solving vertically averaged (or height- 
integrated) radial equations of motion. Thus, steady slim disks neglect the vertical 
structure of flow, and describe essentially flat fluid conflgurations. Although an ex- 
pression for the radial dependence of disk thickness is obtained, the slope of the disk 
surface is usually neglected in the discussion of emergent spectra, where a plane-parallel 
atmosphere is typically considered for the purposes of computing radiative transfer in 
the vertical direction. 

In the case of slim disks, the radial equations of structure are a set of ordinary 
differential equations (Section 4.1). The set of equations is closed by including relations 
describing vertical hydrostatic equilibrium and vertical transport of energy. By solving 
these equations with appropriate boundary (or regularity) conditions one obtains the 
radial proflles of the central temperature, Tc{r), surface density, S(r), height-integrated 
pressure, P{r), half-thickness of the disk, h{r), a radial velocity V{r), emerging flux of 
radiation, F{r), and certain other physical quantities. 

In this Chapter we are taking a flrst step towards constructing truly three-dimensional 
slim disk models, by solving a set of differential equations describing the vertical struc- 
ture of the disk. The resulting 2;-dependence of physical quantities is used to compute 
certain coefficients that enter the radial equations, and that up till now have been 
estimated with algebraic expressions (Eqs. 4.20 - 4.23). We will refer to the resulting 
models as "two-dimensional" (2-D), although it has to be understood that in contrast to 
the thin-disk solutions of (Urpin, 1984; Kluzniak & KiTA, 2000), the actual merid- 
ional flow has not been computed in this paper — the average radial velocity alone has 
been considered in the structure equations. However, the models are now self-consistent 
in the sense that the vertical averages of physical quantities that form the coefficients 
of the radial ODEs do correspond to the vertical structure considered in the radiative 
transfer calculation — in Chapter 4 the vertical structure of the radiative atmosphere 
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was considered a posteriori, and it had no influence on the radial structure of the disk. 



6.1 Radiative transfer 

Accretion disks are not simple rotating fluid configurations which may be described 
by the hydrostatic balance only. The turbulent motions, magnetic fields, mean radial 
velocity, shear and radiation make the description of the disk interior very complicated. 
The simplest approach, which has been successfully used for a long time, is to consider 
disk annuli separately, neglecting the radial gradients, radial net flux of radiation and 
accounting only for the mean flow. Under these conditions, disk layers may be described 
as plane-parallel and are subject to the vertical hydrostatic equilibrium and radiative 
transfer along that direction only. Thus, the problem is reduced to one- dimensional 
and may be relatively easy solved. 

The equation of hydrostatic equilibrium reflects the balance between the gradients 
of pressure and gravitational potential. For accretion disks annuli and cylindrical coor- 
dinates (r, 0, z) it may be put in the form (Eq. 3.12), 

Idp GM 

-— = —z, (6.1) 

p az 

where p = Pgas + Prad- 

The equations describing the radiative part are more complicated. In the follow- 
ing sections we briefly summarize the theory of radiative transfer and introduce the 
Rosseland approximation which will be used in the model described in this Chapter. 

6.1.1 Transfer equation 

The general radiative transfer equation (Rybicki & LiGHTMAN, 1979), 

^ = -a^/j.+>, (6.2) 
ds 

describes the variation of speciflc intensity along the distance s. and stand 
for the absorption and emission coefficients, respectively. Eq. 6.2 takes a particularly 
simple form if another variable, - the optical depth, is introduced. It is defined by 

dr^ = a^ds (6.3) 

and 

ru{s) = [ a,{s)ds'. (6.4) 



If integrated throughout the medium satisfies the condition t,^ > 1, we say that 
medium is optically thick. On the contrary, when < 1, the medium is optically thin 
— the typical photon of frequency u can traverse an optically thin medium without 
being absorbed, while it is not possible in an optically thick medium. 
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Introducing the optical depth, the radiative transfer equation simphfies to 

^ = + S,, (6.5) 

where the source function S^, defined as the ratio of the emission and absorption coef- 
ficients, 

S. = ^ (6.6) 

, has been introduced. 

In general, the intrinsic emission comes from two processes: the reemission after 
absorption (thermal emission) and scatterings. If the former dominates, according to 
the Kirchoff's Law for thermal emission (Rybicki & Lightman, 1979), we have, 

> = aMT) (6.7) 

and 

S, = B,{T), (6.8) 

where By{T) is the Planck function describing the black body radiation (Eq. 4.47). The 
transfer equation is, then, 

^ = _a,(4 - B,{T)) (6.9) 

or 

^ = -/, + i?,(T). (6.10) 
If the medium is dominated by scatterings, we have 

ju = CTuJu, (6.11) 

where is the absorption coefficient of the scattering process (the scattering coefficient) 
and Jj, is the mean intensity of radiation: 

Ju = ^j hdn. (6.12) 

It is then straightforward to get the following equality, 

S, = J, (6.13) 
and the following form of the radiative transfer equation, 

^ = -a,(/,- J,). (6.14) 
as 

If we take into account combined scattering and absorption we have the following 
general form of the transfer equation, 

^ = -a^{I^ - B^{T)) - a^{I^ - J^) = 
as 

= -(a^ + (T^)(/^ - 5^), (6.15) 
with the source function given by, 

OtyB^; -\- (JyJi; , , 

by = . (6.16) 



110 



Vertical structure of slim accretion disks 



6.1.2 Optical depths 

The net absorption coefficient is now equal to + and is called the total extinction 
coefficient. Thus, we should define the optical depth as, 

dry = {a^, + ay)ds. (6-17) 

The mean free path of a photon before scattering or absorption is 

= ^ . (6.18) 

Following the random walk arguments (Rybicki & LiGHTMAN, 1979) we also may 
obtain an expression for the net displacement between the points of creation and de- 
struction of a typical photon. 




(6.19) 



which is also called the effective mean path or the thermalization length. 

From Eqs. 6.18 and 6.19 it follows that the mean path of a photon before scattering 
or absorption may be significantly shorter than the effective mean path. This is the 
case for scattering-dominated media where absorptions are less likely than scatterings. 
The regular total optical depth defined in Eq. 6.17 refiects only the more important 
process. It is therefore convenient to introduce, basing on Eq. 6.19, the effective optical 

depth, 

r* ^ Vra(r, + r,), (6.20) 

where and stand for the optical depths for absorption and scattering, respectively. 

The medium is effectively optically thin, ^ 1, when the effective mean path is 
large compared with the size of the medium — most photons escape out by random 
walking before being absorbed. On the other hand, in effectively thick medium most 
photons thermally emitted at large depths will be absorbed before they get out. The 
medium at large effective optical depths is close to the thermal equilibrium and we have 
both Su By and h ^ (e.g., MlHALAS (1982)). 



6.1.3 Rosseland approximation 

From the previous section we know that deep inside an effectively thick medium the 
radiation is in thermal equilibrium. Let us assume that the properties of the medium 
depend only on the geometrical depth. Symmetry requires that the intensity of radiation 
may depend only on the angle 9 between the given direction and the distinguished 
coordinate. Introducing ii = cos 9, so that 

ds = ^„ = 'A (6.21) 

cosy /i 

we may put the transfer equation (Eq. 6.15) in the following form, 

Iy{z,^,) = Sy-^^^. (6.22) 
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Since at large effective optical depths we have ly = Sy = B^, we may approximate this 
formula with 

Iy{z,fi)^By (6.23) 

The vertical flux of radiation is given by the following integral, 

Fy{z) = j Iy{z,fi)fidQ = 27i j Iy{z,fi)fidfi, (6.24) 

which may be easily calculated. Since only the angle-dependent part contributes, we 
have 

3{a, + ay) or dz- ^^'^^^ 
The total flux may be obtained by integration over all frequencies: 

n.) ^ r FM,. .j-ifr ^^d. (C.20) 

Jo 3 dz Jo ay + ay dT 

Taking into account the following equality, 

dB (T) AaT^ 

^-^du = ^ (6.27) 
and introducing the Rosseland mean absorption coefficient ur. 



1 



(6.28) 



we finally have: 

San oz 

This is the Rosseland or diffusive approximation for the energy flux. It shows that 
the radiative energy transport deep in effectively thick media is of the same nature 
as heat diffusion. It is valid not only in plane-parallel media but everywhere where 
the medium properties change slowly on the scale of the radiation mean free paths. 
The Rosseland approximation says that the vector flux of energy is opposite to the 
temperature gradient and gives its magnitude through Eq. 6.29. 

Heat in gases may be transported not only through radiative processes, but also 
through convection. The latter occurs due to temperature differences which affect the 
density, and thus relative buoyancy, of the fluid. Heavier (more dense) components will 
fall while lighter (less dense) components rise, leading to bulk fluid movement and, as a 
result, transport of heat. The onset of convection may be determined by comparing the 
local values of the thermodynamical radiative and convective gradients. In the model 
described in this Chapter we account for convection by applying the mixing length 
approximation (Paczynski, 1969). 
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6.1.4 Dissipation 

Accretion disk annuli are not conservative in the sense that the flux of energy changes 
with the vertical coordinate due to viscous dissipative processes. Basing on Eq. 3.15 
we expect the whole annuli to generate the following flux of energy at its surface, 

r+h 

F = T{h) = / f^ar^dz, (6.30) 
Jo 

where t^'^ and ar4, are components of the stress and shear tensors, respectively. Using the 
a approach and assuming the heat is not advected, we may write in the non-relativistic 
limit 

F = --aPr^ (6.31) 
2 dr 

with P being the vertically integrated pressure. The models based on vertically averaged 
quantities (like Shakura & Sunyaev (1973) or slim disks) are able to provide the 
total emission but say nothing about the internal distribution of dissipation. So far, no 
analytical model has been established and one has to make ad hoc assumptions or base 
on the results of numerical simulations (e.g., HiROSE ET AL. (2009a)). However, it 
has been shown (S4DOWSKI ET AL., 2009) that the profile of the internal dissipation 
hardly affects the disk observables. 

The simplest and physically justified choice for the dissipation profile is to assume 
that the a prescription is valid throughout the disk interior. Thus, we may take 

^ = ^ap^K, (6.32) 
az 2 

where p is now the local pressure. In the following section we will use the relativistic 
version of this formula. 

However, most recent MHD shearing-box simulations show that the dissipation in 
accretion disks may not follow the local ap dissipation. Surprisingly, the resulting 
maximum of the dissipation rate is not even at the equatorial plane (where both the 
density and the total pressure have maxima). Fig. 6.1 shows the dissipation rate in 
one of the simulations performed by HiROSE ET AL. (2009a). The solid black line 
presents the dissipation rate giving rise to the radiative flux. Its integral is equal to the 
total emitted radiative flux. It is clear that the dissipation profile has two maxima at 
z = ±0.5/i, where h is the disk scale height. As the total (gas, radiation and magnetic) 
pressure peaks at the equatorial plane, the regular ap prescription cannot describe such 
dissipation profile. Nevertheless, as has already been discussed in Section 2.2.2, the a 
formalism works for the total, vertically integrated, quantities. 

6.2 Model 

6.2.1 Basic assumptions, parameters, and coefficients 

We make the same assumptions as for the one-dimensional model (Section 4.1), namely: 
(i) axisymmetry, (ii) stationarity, (iii) Kerr metric, with fixed values of the BH mass. 
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Figure 6.1: Horizontally and time-averaged total dissipation rate (black) and stress 
times shear (dotted) in simulation 1112a of HiROSE ET AL. (2009a). The total dissipa- 
tion rate is the sum of magnetic energy dissipation (red) and kinetic energy dissipation 
(blue). The difference between the stress times shear product and the total dissipation 
rate reflects the fact that some part of the generated heat goes into the mechanical 
energy of gas motion instead of increasing the radiative flux. Figure after HiROSE ET 
AL. (2009a). 
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M, and spin, a, parameters, (iv) symmetry under reflection in the equatorial plane, 
(v) steady mass supply rate, M, through a boundary "at infinity", (vi) zero torque at 
the BH horizon, and (vii) we neglect the loss of angular momentum to both wind and 
radiation, self-irradiation and magnetic pressure. 

A fraction, 1/(1 + /^*^^('")), of the entropy generated locally by dissipative processes 
is released into the radiation field, while the remainder is advected by the gas. 

A unique solution to the slim disk model can only be found if certain additional 
assumptions are made. We make the following arbitrary choice. We neglect the vertical 
variation {z dependence) of the velocity field, considering only its height-averaged value; 
thus, the velocity is always directed radially inwards and is a function of the radial 
coordinate alone. Similarly, we assume there is no z variation of the advection factor 
jadv^^y Dissipation and angular momentum transport are given by the a prescription 
(Shakura & SuNYAEV, 1973), with a constant value of a. We assume that the local 
dissipation rate is proportional to the total pressure, p. For a more detailed statement, 
see Eq. 6.34 and the comment following it. Calculations are carried out for the value 
a = 0.01. 

We are looking for 2-D slim disk solutions at a definite value of mass accretion rate 
for a given Kerr metric. Thus, for a fixed value of a, there are three fundamental 
parameters describing a given slim disk solution: M, a*, and M. 



6.2.2 Vertical structure equations 

We describe the vertical structure of an accretion disk in the optically thick regime by 
the following equations: 

(i) Hydrostatic equilibrium (Kato ET AL., 2008), 

1^ = -nlz, (6.33) 
paz 

where Q± is defined in Eq. 4.13. Although other expressions for the right-hand side of 
Eq. 6.33 can be found in the literature (e.g., Abramowicz ET AL. (1997) includes 
^ 0), the form above is appropriate for our scheme, in which the vertical structure 
is precalculated before any information about the radial variables becomes available. 

(ii) The energy generation equation. We assume that the vertical fiux of energy 
inside the disk J-" is generated according to 

dJ^ 3V f ap \ /A'/V^^ , .^ 

d^= 2c ArTT^J V-^J • ^ ^ ^ 

Strictly speaking, this does not correspond to a constant a prescription, as the term 
(3/2) (M/r^)^/^, which is derived from Keplerian strain, departs somewhat from the 
value that would follow from the actually computed distribution of angular momentum 
(Fig. 6.2). However, as the departure for sub-Eddington accretion rates is small, we 
expect Eq. (6.34) to afford a good approximation. Note that f^'^"' = corresponds to 
the NT disk (going over into the Shakura-Sunyaev disk in the nonrelativistic limit of 
thin disks), /^'^^ > 1 characterizes advection-dominated disks, while f^'^^ < describes 
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those regions where the advected heat is being released. The amount of heat advected 

(iii) Energy transport. The structure of the disk has to be such that the actual 
value of the divergence of the flux corresponds to Eq. (6.34). Radiative transport is 
computed in the diffusive approximation (Eq. 6.29), 

leoT^dT 

^^'^ = dJ' ^'-^'^ 

while convective transport is computed in the mixing-length approximation. Energy 
is transported in the vertical direction through diffusion of radiation or convection 
according to the value of the thermodynamical gradient, which can be either radiative 
or convective. Accordingly, we take 

dlnT ^ r Vrad, for V^ad < Vad /g 
dlnp \ Vconv, for Vrad > Vad 

with the adiabatic gradient given by a derivative at constant entropy: Vad = {d\nT/d\np)s. 
The radiative gradient Vrad is calculated in the diffusive approximation. 



3pkrJ^ 



rad 



(6.37) 



where k/j is the Rosseland mean opacity, and a is the Stefan-Boltzmann constant. 

When the temperature gradient exceeds the value of the adiabatic gradient, we 
have to consider the convective energy flux. The convective gradient Vconv is calculated 
using the mixing length theory introduced by Paczynski (1969) . We take the following 
mixing length, 

H^i = l.OHp, (6.38) 
with pressure scale height Hp defined as (Hameury ET AL., 1998) 



p_ 



Hv = „o2 . • (6-39) 



The convective gradient is defined by the formula 

Vconv = Vad + (Vrad " ^ + w) (6.40) 

where y is the solution of the equation 

7T^Vl/' + ^2/' + w^y-w = Q, (6.41) 
with the typical optical depth for convection = pn^H^x, and w given by 

(Vrad - Vad) • (6.42) 

P 



w'^ V Srrrri / bUa'^T^p VainT 
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The thermodynamical quantities Cp, Vad and [din p/dln.T)p are calculated using stan- 
dard formulae (e.g., Chandrasekhar (1967)) assuming solar abundances {X = 0.70, 
Y = 0.28) and, when necessary, taking the effect of partial ionization of gas on the gas 
mean molecular weight into account. 

We use Rosseland mean opacities (including the processes of absorption and 
scattering) taken from Alexander et AL. (1983) and Seaton et AL. (1994). Fol- 
lowing other authors (e.g., IDAN ET AL. (2008)), we neglect expansion opacities, in 
agreement with our neglect of vertical velocity gradients. 

(iv) We set the following boundary conditions. At the equatorial plane {z = 0) we 
set J^(0) = 0, in accordance with the assumption of the reflection symmetry, while at 
the disk surface, T{h) = 0, we follow the Eddington approximation (MlHALAS, 1982) 
and require J-'{h) = o"Tjj = 2aT^{h). 

In practice, for a fixed r, and prescribed values of Tq and /^*^^, a trial value of 
the central density, pc*, is assumed and the equations are integrated in z until p{z^) = 
10~^^gcm~'^ (as a stand-in for the disk surface, z^, = h). If J^{h) and T{h) fail to 
satisfy the surface boundary condition, the assumed value of pc is adjusted, and the 
integration is repeated, until the condition J^{h) = 2aT'^{h) is met. Convergence is 
usually attained in a few iterations. The emergent flux of radiation at any given r is 
then F = J^{h). 

6.2.3 Radial structure equations 

The radial sector of the model is described by four laws of conservation (of mass, angular 
and radial momenta, and energy), and a regularity condition at the sonic point. These 
equations have exactly the same form as in Section 4.1 with exception to the amount 
of heat advected Q^'^^ which is now calculated precisely basing on the solutions of the 
vertical structure (see the previous section). 

6.3 Numerical methods 
6.3.1 Vertical structure 

The set of ordinary differential equations describing the vertical structure, i.e., Eqs. 6.33, 
6.34 and 6.36 together with appropriate boundary conditions are solved for a given BH 
spin on a three-dimensional grid spanned by the radius r, the central temperature 
Tc, and the advection factor J'^'^^. For a given set of these parameters, we start the 
integration from the equatorial plane [z = 0), and the solution satisfying the outer 
boundary condition is found as described at the end of Section 6.2.2. The resulting 
quantities describing the vertical structure (Tc, S, Q^'^^, P, rji, 772, 773, and 774), together 
with r, are printed out to tables for subsequent use in interpolation routines. As it turns 
out, the first two of these parameters, and E, can be used to uniquely determine all 
the other quantities characterizing the vertical structure, including /^'^^ (see Fig. 6.3). 

Calculating the full grid of vertical solutions for a single value of BH spin takes 
about 5 hours on a 4-CPU workstation. 
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6.3.2 Radial structure 

The radial sector is described by exactly the same differential equations as given in 
Section 4.2. However, the coefficients rji to 774 and the amount of advected heat Q^'^^ 
are now taken from the precise solutions of the vertical structure. 

Similarly like for the one-dimensional case, to obtain a solution one has to solve this 
system of two ordinary differential equations, together with the following regularity 
conditions at the sonic radius r^on, defined by the same conditions as before: 

^f\r =T^\r =0, (6.43) 

' ' son ' ' son ^ ' 

and with outer boundary conditions given at some large radius Tout- The solution 
between the outer boundary and the sonic point is found using the relaxation technique 
(Press, 2002), with treated as the eigenvalue of the problem. The sonic point for 
a^, = is located at 5.9M for accretion rate O.OlMEdd and at 5.0M for 2.0MEdd- We use 
25 mesh points spaced logarithmically in the radius on the section between the sonic 
point and rout = lOOOM. This particular number of grid points is enough to resolve all 
disk features. We have verified that the results are accurately reproduced with a denser 
grid. 

The parameters linked to the vertical structure (Q^*^^, Vi-, V2, Vs^ Vi) given 
E and Tc are linearly interpolated from precalculated tables of the vertical structure 
solutions (for given value of V , S is determined directly from the mass conservation). 
The radial derivatives dln?7i/dlnr, dln?72/dlnr, dln?73/dlnr, and dln?74/dlnr are 
evaluated numerically from the ?7i, 772, r/s, and 7/4 profiles in the previous iteration step. 
A relaxed solution is obtained in a few iteration steps. Once a solution outside the sonic 
point is found we numerically estimate the radial derivatives of V and Tq at the sonic 
point using values given at r > r^^n and use these derivatives to start direct integration 
inside the sonic point. 

The solution thus obtained may then be used as a trial solution when looking for 
the relaxation solution of another slim disk, i.e., when one of the three fundamen- 
tal parameters (M, a,,, M) has a slightly different value. Each relaxation step takes 
approximately 5 seconds on a single-CPU workstation. 



6.4 Disk structure 

In the following two subsections we present and discuss both the radial and vertical 
structure of the two-dimensional model of slim accretion disks. All the solutions, if not 
stated otherwise, were computed assuming a = 0.01 and M = 10 M©. 

6.4.1 Disk radial structure 

The radial profiles are very similar to those described in Section 4.3 (detailed comparison 
is given in Section 6.5). However, to make this section complete, we shortly discuss 
them again. In addition to figures presented before, we show disk thickness and surface 
density dependence on BH spin (Fig. 6.4) as well as profiles of the total and effective 
optical depths (Figs. 6.6 and 6.7). 
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Figure 6.2: Angular momentum (u^) in a Schwarzschild slim disk for two accretion 
rates. For very low accretion rates the angular momentum follows the Keplerian profile 
(dotted line) down to the ISCO. For high accretion rates the flow is super-Keplerian 
between the "center" of the disk at Vccn and the "potential spout" at rpot- The vertical 
dot-dashed line on this and subsequent figures denotes the location of the ISCO. 
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Angular momentum 

The angular momentum profiles for our accretion disk solutions near a non-rotating BH 
are presented in Fig. 6.2. Results for two values of mass accretion rate are shown, M = 
O.lMEdd and 2.0MEdd- For the lowest accretion rates the profiles follow the Keplerian 
profile and reach its minimal value {Cin in Eq. 4.10) at the marginally stable orbit 
(ISCO). The higher the accretion rate, the stronger the departure from the Keplerian 
profile. For low values of a (e.g., 0.01) the disk is sub-Keplerian at large distances and 
super-Keplerian at moderate radii. The Keplerian profile is crossed again at a point 
located inside the marginally stable orbit, and corresponding to what is usually called 
"the cusp" or "the potential spout". For a detailed study of the physics of differently 
defined inner edges of accretion disks see Section 4.4. 



S-curves 

Figure 6.3 presents slim disk solutions at r = 20M on the Tc-T, plane, for a non-rotating 
BH. Solutions of the polytropic, height-averaged models are presented for comparison 
(for detailed discussion see Section 6.5). The locus of solutions for various values of the 
mass accretion rate has the shape of the so-called "S-curve" (Abramowicz et AL., 
1988). The lower, gas-pressure dominated branch accurately follows the track of radia- 
tively efficient solutions {f'^'^^ = 0). The middle, radiation-pressure dominated branch 
is reached at M ^ O.lMsdd- As advection becomes significant, the slim disk solution 
leaves the /^"^^ = track and moves to higher advection rates. Around M = 5MEdd, 
the solutions enter the upper advection-dominated branch corresponding to /'^'^^ > 1.0 
(more than 50% of heat stored in the accreted gas). At M = 20MEdd this rate increases 
almost up to 80% {f^'^^ ^ 4.0). 



Surface density 

Profiles of the surface density for 2-D slim disk solutions for a non-rotating BH are 
presented in the upper panel of Fig. 6.4. Different regimes, corresponding to different 
branches of the "S-curve" on the (S, Tq) plane are visible. For large radii the sur- 
face density increases with increasing accretion rate (the lower gas-pressure dominated 
branch), while this relation is opposite for moderate radii (the middle radiation-pressure 
dominated branch). For accretion rates M > 5.0MEdd the upper advection-dominated 
branch would be reached. The local maxima in the surface density profiles (discussed in 
detail in, e.g., S4DOWSKI (2009)) are visible for moderate accretion rates (~ 0.5MEdd)- 
The bottom panel of Fig. 6.4 presents corresponding profiles of the radial velocity V as 
measured by an observer co-rotating with the fluid. 

The surface density dependence on BH rotation is presented in Fig. 6.5. The profiles 
are shifted to lower radii as the inner edge of the disk moves inward for higher BH spins. 
The outer parts of the accretion disk are insensitive to the metric. 
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Figure 6.3: The Tc-S plane at r = 20M for a non-rotating BH (a^, = 0). The dotted 
hnes connect solutions for the vertical structure of slim disks that have the same value 
of the advection parameter J^'i^. The locus of standard (radiatively efficient, /^'^^ = 0) 
disk solutions is shown with the thick dotted line. The solid thick line represents the 
vertical slim disk solutions for different accretion rates (indicated by triangles), and 
the dashed line presents corresponding solutions of the conventional polytropic slim 
disk model (see Section 6.5). The difference between the two lines in the low M limit 
corresponds to the difference in S between the two models (see Fig. 6.18). 
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Figure 6.4: Profiles of tlie surface density (upper panel) and corresponding values of 
radial velocity V (bottom panel) of a slim disk for a non-rotating BH. Solutions for 
different accretion rates are presented. 
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Figure 6.5: Profiles of tlie surface density for slim disks at a constant accretion rate 
(M = O.lMEdd) and various BH spins. 



Optical depth 

In the top panel of Fig. 6.6 we plot the total optical depth of the vertical slim disk 
solutions for different accretion rates and a* = 0. The total optical depth, 

Ttot = / i^Rp dz, (6.44) 
Jo 

where the total opacity coefficient kr, which includes the processes of absorption and 
scattering, is closely related to the surface density. Indeed, the radial profiles of the 
optical depth shown in Fig. 6.6 follow the corresponding profiles of surface density. Any 
differences in the profiles come from the dependence of the opacity coefficient on local 
density and temperature. Outside the ISCO the total optical depth is always large 

(Ttot > 10^). 

The diffusive approximation for radiative transport may only be used if photons are 
absorbed, otherwise LTE cannot be established. In a scattering-dominated atmosphere, 
the effective optical depth is then the relevant quantity to be used in checking for the 
self-consistency of the diffusive approximation. The bottom panel of Fig. 6.6 presents 
corresponding profiles of the effective optical depth, which is estimated in the following 
way: 

rh 

Tcff = / ^y{nR - Kes)l^R P dz (6.45) 
^0 

where k^s = 0.34 cm^g"^ is the mean opacity for Thomson electron scattering. For 
M > O.lMEdd the inner region of the disk becomes radiation-pressure dominated (it en- 
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Figure 6.6: Optical depth for a = 0.01 slim disks around a Schwarzschild black hole at 
different accretion rates. Top panel: The total optical depth as a function of the radius. 
Bottom panel: The effective optical depth as a function of the radius. The ISCO is 
shown at r = 6M. The shaded region in the bottom plot indicates the region where 
the diffusive approximation is invalid. 
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Figure 6.7: Profiles of tlie effective optical depth of a Schwarzschild slim disk for three 
values of viscosity (a = 0.01, 0.05, and 0.1), calculated for two accretion rates, O.lMEdd 
(dashed lines) and l.OMEdd (solid lines). 

ters the middle branch on the corresponding S-curve), the surface density decreases with 
increasing accretion rate, and electron scattering begins to dominate absorption. There- 
fore, the effective optical depth decreases with increasing accretion rate and reaches 
values Tcfj < 1 in the inner parts of the disk for accretion rates above l.OMEdd- As a 
result, for M > l.OMEdd? the diffusive approximation can no longer be applied, and 
our present approach to solving for the disk structure breaks down. In Fig. 6.7 we 
exhibit the dependence of the effective optical depth on the value of a. In general, 
Teff is inversely proportional to a (as is the surface density). At lower accretion rates 
(M < O.lMEdd) the effective optical depth remains large even for high values of a. 

Flux profiles 

Profiles of the flux emitted from the disk surface {F = o"Tjj) in the case of a non-rotating 
BH are presented in Fig. 6.8. Results corresponding to accretion rates from 0.01 up to 
S.OMEdd are shown. For the lowest rates the emission from inside the marginally stable 
orbit is negligible as expected in the standard accretion disk models. This is no longer 
true for higher accretion rates, and the advection of energy causes significant emission 
from smaller radii. For super-Eddington accretion rates the emitted flux continues to 
grow with a decreasing radius even inside the marginally stable orbit. Radiation coming 
from the direct vicinity of the BH is suppressed by the gravitational redshift (the g- 
factor, see Section 4.5.2). Therefore, an observer at infinity will observe a maximum in 
the profile of the effective temperature even for the highest accretion rates. 
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Figure 6.8: Flux emitted from the surface of a slim disk at five accretion rates onto a 
Schwarzschild black hole. At high accretion rates significant emission from within the 
ISCO is clearly visible in the figure. 

The increase of the advective flux with increasing accretion rate is clearly visible in 
Fig. 6.9. The ratio of the heat advected to the amount of energy emitted is presented 
for different accretion rates. For very low accretion rates these profiles approach the 
limit of a radiatively efficient disk {f^'^^ = 0), and the advection component becomes 
significant for higher accretion rates. Some part (up to 30% at r = 20M for 2.0MEdd) 
of the energy generated at moderate radii is advected with gas and radiated away at 
r < lOM. This causes the significant change in the emitted flux profile at the higher 
accretion rates visible in Fig. 6.8. 

In Fig. 6.10 we present the emitted flux profiles for different BH angular momenta at 
a constant accretion rate M = O.lMEdd- These profiles coincide at large radii where the 
influence of the BH rotation is negligible; however, the higher the BH spin, the closer 
to the horizon the marginally stable orbit. Therefore, in the case of rotating BHs, 
the accreting matter can move much deeper into the gravity well, compared to non- 
rotating BHs. This effect leads to an increase in the disk luminosity and hardening of its 
spectrum, which can be inferred from Fig. 6.10 — the higher the spin, the higher the disk 
luminosity, and the higher the flux (which corresponds to the effective temperature). 

Photosphere location 

The flux observed at inflnity may be obtained by performing ray tracing of photons 
emitted from the accretion disk (see Section 6.5). Scattering in the layer above the 
photosphere must also be taken into account. An accurate calculation requires de- 
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Figure 6.9: Profiles of tlie advection coefficient / for different accretion rates 
(Scliwarzscliild black hole). 




Figure 6.10: Flux profiles at fixed accretion rate [M = O.lMEdd) for five values of BH 
spin. 
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Figure 6.12: Profiles of pfiotospheric height at a constant accretion rate (M = O.lMEdd) 
and various BH spins. 
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tailed knowledge of atmospheric properties, including the location of the photosphere. 
This is particularly important when accretion rates are high and the disk is no longer 
geometrically thin (S4DOWSKI ET AL., 2009). 

In Fig. 6.11 we plot the profiles for the 2;- location of the photosphere, /iphot? ob- 
tained in our model at different accretion rates for a non-rotating BH. Clearly, for 
high accretion rates, M > O.lMEdd, the inner regions become thicker (effects of radi- 
ation pressure). For M = l.OMEdd, the ratio /iphot/''" reaches values as high as 0.25. 
Close to the ISCO the height of the photosphere rapidly decreases because of vig- 
orous cooling. Although the rapid change in disk thickness violates the assumption 
of the hydrostatic equilibrium that we make when solving for the disk vertical struc- 
ture, the accelerations connected with the vertical motions involved are much lower 
than the vertical component of gravity. Indeed, the vertical accelerations are close to 
Vrdvz/dr ~ Vrd{vrdh/dr)/dr ~ Vgd{h/r)/dr ~ rfl^lh/r)^ ~ Q\h{h/r). In deriving this 
estimate we liberally assumed that dh/dr ~ 1 and used the fact that the rapid decrease 
in disk height occurs near the sonic point, while the speed of sound Vg is approximately 
rQ±{h/r). Thus the acceleration terms modifying Eq. 6.33 would be smaller than the 
gravitational acceleration by a factor of a few percent: [h/r) ~ 10~^. In Fig. 6.13 we 
plot the dynamical and gravitational components of the vertical equilibrium equation 
at the photosphere. It is clear that the former is at least 10 times smaller than the lat- 
ter at the sonic radii for M < MEdd- Therefore, in all likelihood, our results correctly 
describe the disk structure that would be obtained without assuming strict hydrostatic 
equilibrium. 

In Fig. 6.12 we present radial photosphere profiles, at a fixed accretion rate and 
different values of the BH spin. The photospheric heights coincide for large radii. In 
the inner regions of the disk, the height of the photosphere increases with BH spin, 
reflecting the increased luminosity and radiation pressure. 

6.4.2 Vertical structure 

In Fig. 6.14 we present a vertical cross-section of a Schwarzschild slim disk for M = 
O.OlMsdd- At this accretion rate the disk is radiatively efficient and no advection of 
entropy is expected. The top panel presents the radial profiles of the photosphere and 
the disk surface (defined as a layer with p = 10~^^g/cm^). 

The total optical depth reaches values as high as ~ 20000 on the equatorial plane and 
decreases monotonically towards the photosphere at r = 2/3 (Fig. 6.14, second panel 
from the top). Within the marginally stable orbit, the total optical depth significantly 
decreases, as shown in Fig. 6.6. 

The third panel of Fig. 6.14 presents the temperature, which decreases with height 
from Tc at the equatorial plane to T^s at the photosphere. The maximum value of the 
temperature is attained on the equatorial plane at r ~ lOM. 

Density is presented in the fourth panel, and the next panel presents the vertical flux 
that is generated inside the disk according to Eq. 6.34. It is set to zero on the equatorial 
plane by reflection symmetry, and then rapidly increases, because in an alpha disk the 
dissipation is proportional to the pressure, which reaches its maximum on the z = 
plane. Close to the disk surface, where pressure is almost negligible, the flux slowly 
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Figure 6.13: Comparison of the dynamical (tvdf^/dr ^ Vd{VdH/dr)/dr) and gravita- 
tional {Q^H) components of the vertical equilibrium equation at the photosphere. The 
stars denote locations of the sonic radii. 

settles down to the emitted value. At the accretion rate chosen for the figure the flux 
rapidly decreases inside the marginally stable orbit. 

The bottom panel of Fig. 6.14 presents the non-monotonic distribution of the ther- 
modynamical gradient (Eq. 6.36), which ranges between 0.2 and 0.4. Therefore, the 
disk's vertical structure cannot be described by a simple polytropic relation. Moreover, 
in the region of the highest temperature (close to the equatorial plane at moderate 
radii, r ~ 20M), the heat is transported upward through convection. 

The vertical structure of an accretion disk with ten times higher accretion rate, 
M = O.lMEdd; is presented in Fig. 6.15. The general picture remains the same, since 
the advection of heat is still insignificant. However, as the temperature increases, the 
inner disk regions become dominated by radiation pressure. For O.lMEdd the convective 
region extends from the marginally stable orbit up to 300M and covers more than half 
of the disk thickness. 

The disk structure is significantly different in the case of a high accretion rate (e.g., 
l.OMEdd); with a significant amount of advection. The vertical cross- sect ions of the slim 
disk are presented in Fig. 6.16. The inner regions are dominated by radiation pressure, 
so the disk geometrically thickens and the photosphere is now higher. The total optical 
depth mostly follows the surface density (compare Fig. 6.6), and therefore decreases 
considerably towards the BH in the inner parts of the disk. 

The temperature maximum is again located at the equatorial plane close to r 
lOM. The maximum of the effective temperature (corresponding to the vertical flux 
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Figure 6.14: Vertical structure of a slim disk for M = O.OlMEdd and a* = 0. The 
top panel presents the surface of the disk (green dashed line), and the photospheric 
surface (red solid line). The other panels present the structure of the disk below the 
photosphere. Top to bottom: total optical thickness, temperature, density, vertical flux 
of energy, and the thermodynamical gradient. The blue solid line in the bottom panel 
delimits the convective region. 




Figure 6.15: Same as Fig. 6.14 but for M = O.lMEdd- 




Figure 6.16: Same as Fig. 6.14 but for M = l.OMEdd- 
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Figure 6.17: Meridional profiles of density for three accretion rates: M = 0.01 (top), 
0.1 (middle) and l.OMEdd (bottom panel) in (r,z) coordinates. The violet boundaries 
show the location of the photosphere. The black hole is described by Mbh = IOMq, 
and a^: = 0. 



shown in the fifth panel (compare also Fig. 6.8) is shifted inwards, down to r ^ 8M. 

The fourth panel of Fig. 6.16 presents the density distribution. Despite the fact that 
the surface density monotonically increases outward (see Fig. 6.4), p has two maxima in 
the equatorial plane: at r 200M and r f« 6M. Finally, the bottom panel presents the 
thermodynamical gradient distribution. For such a high accretion rate, the convective 
zone extends nearly to the photosphere for r < 200M and is present up to r = 2000M. 

In Fig. 6.17 we present the density in the meridional plane for three different ac- 
cretion rates. The equatorial plane lies in the middle of each plot. The violet bound- 
aries denote the photosphere. The two maxima of density at the equatorial plane for 
M = l.OMEdd are clearly visible also in this representation, as are some other features 
discussed above. 



6.5 Comparison with height- averaged shm disk solu- 
tions 

In this section we compare the two-dimensional (2-D) slim disk solutions, where the 
radial structure equations are coupled to those for the vertical structure, with the 
standard polytropic slim disk model, introduced in Chapter 4, in which the slim disk 
equations and properties are averaged over the thickness of the disk. The latter may 
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be extended to mimic the more detailed 2-D solutions by introducing two additional 
parameters: fp and fn- 

Under the one- zone approximation and assuming that radiation is transported in 
the vertical direction only through diffusion, one obtains the following formula for the 
total flux emitted from disk surface (compare Eq. 4.17), 

= /i^^, (6.46) 

where factor fp has been introduced to account for inaccuracies arising from this ap- 
proximation, as well as from the dominance of the disk convection in certain regions 
(as discussed in §6.4.2). Now, advective cooling takes the form 

. .^P^f _ f/^. (6.4T) 

r-^ ar Shn 

One has to remember that disk thickness, h, introduced in Eq. 4.12, is not the exact 
location of the photosphere but rather disk pressure scale. Therefore, we introduce a 
factor fh relating these quantities, 

Vot = fhh. (6.48) 

Usually (like in Chapter 4), the following values of N (the polytropic index), fn 
and fp are assumed: 

N = 3.0, 

fH = 1.0, (6.49) 
fp = 1.0. 

In Fig. 6.18 we compare radial profiles of the flux, photospheric height and surface 
density of the 2-D slim disk model described in this Chapter (consistently taking its ver- 
tical structure into account) with profiles obtained for the polytropic, height-averaged, 
slim disk, as presented in Chapter 4. The comparison was carried out for two accre- 
tion rates (0.1 and l.OMEdd) and a = 0.01. For the lower accretion rate (O.lMgdd)) 
the disk is radiatively efficient (advective cooling is negligible), and therefore both flux 
proflles almost coincide and correspond to the NOVIKOV & Thorne (1973) solutions. 
However, the disk photosphere location in the polytropic model turns out to be overes- 
timated by more than 20% in the region corresponding to the maximal emission. The 
proflles of the surface density do not coincide either, the 2-D solutions giving values 
~ 25% lower than the corresponding polytropic solutions (compare Fig. 6.3). 

For the higher accretion rate (LOMedd)? the flux proflles remain similar (up to 1%). 
However, as advection becomes important, the emission is shifted inwards with respect 
to the Novikov & Thorne proflle. The photosphere location in the polytropic model is 
overestimated by ~ 30% and the surface density by ~ 20%. 

A question arises as to whether such differences in the flux, photosphere, and sur- 
face density proflles affect the resulting disk spectrum. In Fig. 6.19 we present spectral 
proflles and their ratios (2-D to polytropic) for two accretion rates. The spectra were cal- 
culated with ray-tracing routines (BURSA, 2006) using the BHSPEC package (Davis 
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Figure 6.18: Comparison of the flux, disk thickness and surface density profiles cal- 
culated using the 2-D (this Chapter, solid lines) and the usual polytropic (Chapter 4, 
dotted lines) slim disk models for a = 0.01. The solutions for two accretion rates (0.1 
and l.OMEdd) are presented with thick and thin lines, respectively. 
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Figure 6.19: The upper panel presents spectral profiles of the 2-D and polytropic solu- 
tions for two accretion rates (0.1 and l.OMEdd), at inclination angle i = 70° and distance 
d = lOkpc. The bottom panel presents ratios of the corresponding spectra (of 2-D to 
polytropic solutions) for both accretion rates. 



ET AL., 2005), assuming the inclination angle i = 70° and distance to the observer 
d = lOkpc. As BHSPEC gives tabulated solutions of the full, frequency-dependent, 
radiative transfer equations for the disk vertical structure taking the Compton scat- 
terings into account in the disk atmosphere, the spectra presented in Fig. 6.19 are 
not those of a simple multicolor blackbody. However, one should be aware that using 
BHSPEC for calculating spectral color correction is not consistent with the assumed 
vertical structure (calculated or height-averaged) because it is based on a stand-alone 
disk model. 

The general shape of the spectra is similar for both types of slim disk models, 
because the emission profiles nearly coincide. However, the spectra are not identical. 
The bottom panel of Fig. 6.19 presents ratios of the spectral profiles of the corresponding 
solutions for each accretion rate. They all coincide att low energies (< O.lkeV), while 
for higher energies the discrepancies are as large as 3% for l.OMEdd at 5keV. These 
differences are attributed to (slightly) different profiles of the flux, the photosphere, and 
the surface density. We conclude that the proper treatment of the vertical structure 
hardly affects the spectra of slim disks for this range of accretion rates (M < Msdd) 
and the viscosity parameter (a < 0.01). 
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We end this section by giving fitting formulae for A^, fn, and fp, which approximate 
the full numerical 2-D solutions described in this paper with a polytropic slim disk 
model. An advantage of using these formulae lies in avoiding the need to perform 
time-consuming calculations of the vertical structure and avoiding numerical problems 
connected to interpolation in the vertical solutions grid. The formulae for the polytropic 
model parameters for a = 0.01 are 



A^ 


= 3.25 X Sn, 




fn 


= 0.63 X Sh, 


(6.50) 


fp 


= 0.94 X Sf, 





where the spin correction coefficients Sn, Sh, and Sp are given by 

= 1 + 0.002 (6 - risco/M), 
Sh = 1 + 0.003 (6 -risco/M), 
Sp = 1 + 0.064 (6 - risco/M). 

Here, risco is the radius of the marginally stable orbit. For a non-rotating BH, at the 
radius r = 7M (corresponding to the highest disk effective temperature at M = MEdd), 
the fitting formulae are accurate to 1% for the emitted fiux, the photospheric height 
and the surface density. 



6.6 Discussion 

Motivated by a desire to explain and fit the observed spectra of accreting BHs in binary 
systems to theoretical models, we have developed a 2-D model of optically thick slim 
disks. These should be particularly relevant to transient binaries. In quiescence, their 
inner disk regions are described well by optically thin advect ion-dominated accretion 
fiows (ADAFs; see e.g., Lasota et AL. (1996); DuBUS ET AL. (2001)). However, a few 
BH systems (e.g., GRS 1915+105 and LMC X-3) have been observed in thermal states 
corresponding to disk luminosities higher than 0.3LEdd (McClintock & Remillard, 
2003; Steiner et AL., 2010), and modeling these require optically thick models going 
beyond the standard thin disks. 

In this Chapter we presented a 2-D slim disk model, in which the radial and vertical 
structures are coupled. Such an approach eliminates arbitrary factors that infiuence 
solutions of the usual polytropic slim disk model. The results were obtained under two 
key assumptions: an alpha disk was assumed (dissipation proportional to pressure), 
with a uniform value of a, and the fraction of the generated entropy that is advected 
was computed at every radius under the assumption that this fraction does not vary with 
the height above the disk plane (Eq. 6.34). Both of these assumptions seem arbitrary, 
and we can offer no physical motivation for the (conventional) choice we made. 

Under these assumptions and for the value a = 0.01 of the viscosity parameter, 
we computed and presented the detailed structure of 2-D slim disks, parametrized 
by the mass accretion rate, and the two Kerr metric parameters, M and a. Somewhat 
surprisingly, the spectra observed at infinity from such disks differ by only a few percent 
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(and only at high energies) from those obtained from previously considered slim disk 
models (in which the equations and structure correspond to a height average over a 
polytropic atmosphere). Such differences are unlikely to introduce any large corrections 
to spin measurements based on X-ray continuum fits made with corresponding height- 
averaged polytropic models of slim disks. Therefore, applying the spectral routines 
based on regular, polytropic slim disk solutions (e.g., slimbb, described in Section 4.5.2) 
is justified. 

One has to be aware that the model of vertical structure presented here is only an 
approximation of the real physical processes taking place in disk interiors. The diffusion 
approximation and the convection treatment in the mixing length approach are known 
to successfully describe media with large effective optical depths but break down when 
the disk becomes optically thin. We have shown that the effective optical depth of slim 
accretion disks may drop below unity for super-Eddington luminosities and sufficiently 
high values of a. For such conditions, a more sophisticated model of radiation transfer 
should be implemented. However, for a < 0.01 and L < L^^dd the assumptions of this 
work are self-consistent. For higher values of a their range of applicability is limited 
to lower luminosities (e.g., to O.SLgdd for « = 0.1). Kerr-metric slim disks with low 
effective optical depth were discussed by BeloborODOV (1998), who finds them to be 
significantly hotter than the optically thick ones. 

Another remark is connected to the fact that one can expect winds to be blown out of 
the disk surface at super-Eddington luminosities. Such a phenomenon may significantly 
change the disk structure, e.g., its thickness. This feature of slim disks, not described 
in our calculations, has recently been recently cleverly modeled by DOTANI & Shaviv 
(2010). 



Chapter 7 

Self-irradiated slim accretion disks 



The slim disk models presented and discussed in Chapters 4, 5 and 6 make a few funda- 
mental assumptions which make the problem of solving the hydrodynamical equations 
simpler, e.g., the negligence of the angular momentum taken away by photons in the an- 
gular momentum equation (Eq. 2.16) makes its direct integration possible and puts the 
hydrodynamics into the form of ordinary differential equations instead of an integral- 
differential problem. The same physical process may be accounted for consistently in 
the standard model of a relativistic thin and radiatively efficient disk (NOVIKOV & 
Thorne, 1973). 

The returning flux of radiation and angular momentum due to the self-irradiation 
in a curved spacetime is another factor which has been so far neglected. Again, only 
in the case of thin disks one can account for this effect in a (relatively) simple way. It 
has been shown (Ll ET AL. (2005), Kato ET AL. (2008)) that for thin disks with no 
torque at their inner edge the returning radiation has little influence on the inner disk 
regions and therefore on the emerging flux and high-energy disk spectrum. However, 
as was shown in the previous Chapters, accretion disks with high accretion rates are 
geometrically thick and the effect of self-irradiation may be much more significant. 
Unfortunately, including the appropriate terms into the disk equations once again leads 
to a complicated set of integral-differential equations. Therefore, a different approach 
for solving them is necessary. 

In this Chapter we present a state-of-art numerical method of solving slim disk 
equations including all of the above-mentioned factors. We limit our discussion to the 
case of a non-rotating BH. More detailed study will be included in a separate paper. 

7.1 Equations 

We start with presenting modified forms of the vertically integrated equations describing 
an advective accretion disk which include the flux of angular momentum taken away 
by photons and the fluxes of radiation and angular momentum coming back with the 
returning radiation. Only two equations need to be modified: the angular momentum 
conservation and the energy balance. Deriving them we assume that the disk albedo 
is zero, i.e., all the returning radiation is absorbed and thermalized by the disk. We 
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note that this assumption may not be satisfied in the disk regions which are effectively 
optically thin (e.g., deep inside marginally stable orbit). Our approach gives the upper 
limit on the impact of the returning radiation on the disk structure and neglects the 
reflected radiation. 



7.1.1 Angular momentum conservation 

Let us take the angular momentum conservation law in the general form (Eq. 2.16) 
assuming stationarity [dt = 0) and leaving the general form of the (r, 0) component of 
the stress-energy tensor, 

MdC d , 

2^ d^ = -d;:^^ ^'-'^ 

This expression neglects the angular momentum fluxes due to the emerging and return- 
ing radiation which are equal, 

Fcout = CF, (7.2) 
Fc,in= I I £(rOF(r')7^(r,r^p)dfidr^ (7.3) 

J disk surface J hemisphere 

respectively. F = F{r) denotes the emerging flux at given point on the disk surface, 
C = C{r) is the corresponding angular momentum, 7^(r, r', p) is a ray-tracing function 
determining if a photon emitted from r' along the direction p returns to the disk at 
radius r, and which accounts for the change of the photon energy (expressed by the 
g-factor, see Section 4.5.2 or Ll ET AL. (2005)). 
Including these terms leads to, 

MdC d , „ , ^ ^ 

7r:r = -^t^^ t,.^) + ^^c^o.t - rFc^^ (7.4) 

Zir ar ar 
and further, after radial integration, to, 

^{C- An) = r\-Tr^ + Tout - Ti„), (7.5) 

ZTT 

where 

Tout = 4 / ^'^/:,outdr' (7.6) 

and 

Tin = 4 / r'Fcindr'. (7.7) 

Calculating these integrals requires a priori knowledge about the emission and angular 
momentum proflles between the horizon and given radius r. We will solve this prob- 
lem using the relaxation method where the integrals are calculated basing on previous 
relaxation steps. Details are given in Section 7.2. 



7.2 Numerical methods 
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7.1.2 Energy balance 

The energy balance equation (Eq. 4.17) relates the amount of the advected heat (left 
hand side) to the difference between the heating and cooling (right hand side). Origi- 
nally, the latter involved the viscous heating and radiative cooling only. When account- 
ing for the returning radiation one has to include an additional heating term as the 
incoming flux of radiation brings its energy and stores it in the disk heating it up. The 
expression for the advected heat takes the following form, 



where 



is the viscous heating. 



gadv _ gvis _ grad ^ grct ^7 g^, 

Q- = -aP^f (7.9) 
dr 

0- = ^ (7.10) 

is the radiative cooling, and the returning flux Q'^^^ is given by, 

Q"''= I I F(r')7^(r,r',p)dfidr'. (7.11) 

J disk surface J hemisphere 

The integral will be calculated using the emission profile from the previous iteration 
step as described in detail in the next section. 



7.2 Numerical methods 

Including self-irradiation into the slim disk model introduces integrals which cannot be 
directly evaluated because they depend, e.g., on the emission and angular momentum 
profiles which are not known at the time of solving the equations. However, the exact 
solution may be found using an iteration scheme where these integrals are calculated 
basing on the profiles from previous iteration steps. The solution is found when and if 
such a scheme converges. Fortunately, this is the case for self-irradiated slim disks. In 
this paragraph we give details of the numerical methods used to solve the problem. 

We base on the equations for the radial structure of the polytropic slim accretion disk 
described in Chapter 4. We take without any modifications the equation of continuity 
(Eq. 4.7), the radial momentum balance (Eq. 4.8), the vertical equilibrium formula 
(Eq. 4.12) and the expression for the advected heat (Eq. 4.16). We modify the formula 
for the angular momentum conservation and the heating-cooling balance and take the 
forms derived in the previous two paragraphs, 

^{C- An) = (-T,^ + e(Tout - Ti,)) , (7.12) 
Ztx 

gadv _ gvis _ grad ^ ^gret (-^^^3) 

introducing a dimensionless parameter ^ satisfying < ^ < 1. 
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In the beginning we put Tout = Tin = Q'^^^ = and solve the set of regular slim disk 
equations. Basing on this solution, we calculate the integrals in Eqs. 7.6 and 7.7 as 
well as the returning radiation profile Q''^* (Eq. 7.11). We set ^ to some small number 
(typically 0.001) and relax the equations (in the same way as described in Section 4.2) 
to obtain a transonic solution. Once the transonic solution is found we slightly (not to 
lose convergence) increase value of ^ making the impact of the self-irradiation effects 
stronger. During the relaxation in ^ we keep the calculated values of Tout, ^in and Q'^'^^ 
fixed. This procedure is repeated until the solution with ^ = 1 is found. Then, new 
profiles of the integrals and returning radiation are calculated. Basing on them, we 
relax the solution from the one obtained in the previous set of iterations to a new one, 
following the most recent estimate of Tout, Tin and Q^^^. This scheme is repeated until 
a converged solution is found. 

Usually, around 50 iterations in C, are required and the whole scheme has to be 
repeated few times. Due to the numerical evaluation of the self-irradiation, the re- 
quired computational time increased significantly (up to 10 minutes on a single-CPU 
workstation) in comparison with the regular, non-irradiated disk model. 

7.3 Results 

In this section we discuss the most important features of self-irradiated solutions of 
slim disks for a non-rotating BH. More detailed study, including discussion of the im- 
pact of self-irradiation on accretion disks around highly-spinning BHs, will be put in a 
forthcoming paper. 

7.3.1 Emission profiles 

The profile of emission, corresponding to the rate of radiative cooling, is the most 
important factor determining the accretion disk spectrum. Other factors which are im- 
portant for the spectrum are: location of the disk photosphere (the location where the 
ray-tracing starts) and disk surface density (related to the intrinsic hardening of the 
spectrum due to photon up-scattering by hot electrons, see Section 4.5.1). These addi- 
tional factors, as our study shows (e.g.. Fig 7.1), are slightly affected by self-irradiation, 
and the changes in the emission profile dominate. Thus, we put most attention to 
discussing the impact of returning radiation on the emergent flux profiles. 

In Fig. 7.2 we plot profiles of the emitted and returning fluxes of radiation for slim 
disks with three accretion rates: 0.1, 1.0 and lO.OMsdd (marked with different colors). 
The dotted lines denote solutions of the regular, non-self-irradiated model presented in 
Chapter 4. The dashed lines give the profiles of the returning radiation, while the solid 
lines show the emission from the new, self-irradiated solutions. 

The disk with the lowest accretion rate presented (O.lMgdd) is very thin {h/r < 0.1, 
see Fig. 7.1) and therefore, the returning radiation is a result of the light bending only — 
there is hardly any direct irradiation by the disk itself and curved geodesies are required 
to direct photons back to the disk. Thus, the returning radiation for this accretion rate 
(blue, dashed line) with good accuracy corresponds to the profile of returning radiation 
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Figure 7.1: Disk thickness profiles of the self-irradiated solutions (solid lines) for three 
accretion rates: 0.1 (blue), 1.0 (orange) and lO.OMsdd (green), a = 0.1 and a non- 
rotating BH. The dashed lines (which significantly differ from the solid lines only for 
the highest accretion rate) present corresponding profiles for the standard slim disk 
model (Fig. 4.6). 
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Figure 7.2: Profiles of the emitted and returning fluxes of radiation (solid and dashed, 
respectively) predicted by the self-irradiated slim disk model for a non-rotating BH and 
a = 0.1. Profiles of the emitted flux for regular (non-irradiated) disks are also presented 
for comparison (dotted line). Colors correspond to different accretion rates: 0.1 (blue), 
1.0 (orange) and lO.OMEdd (green). 
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calculated for the standard thin disk (Ll ET AL., 2005). For the case of a non-rotating 
BH, the profile of returning radiation is roughly a power law with the same exponent 
as the emission at r ^ HscO) i-C-? -^ret oc r~^. At large radii the returning radiation 
is equal to roughly 3% of the radiation emitted. The relation between the emitted 
and returning fluxes changes rapidly inside the marginally stable orbit. Inside this 
particular radius, there is little viscous heating, and therefore the standard solutions 
emit hardly any radiation. On the other hand, the radiation emitted from outside 
the marginally stable orbit is efficiently focused by the BH on the regions close to the 
horizon. Therefore, the returning flux increases with decreasing radius, and dominates 
over the viscosity-generated emission. 

The solid blue line presents the ultimate solution for O.lMEdd, which takes into 
account the impact of self-irradiation and the angular momentum taken away by pho- 
tons, and which has relaxed following the numerical procedure described in Section 7.2. 
Outside the marginally stable orbit the disk is heated up by viscosity, while the self- 
irradiation dominates inside. One can expect that an accretion disk at such low accre- 
tion rate is radiatively efficient. This is in fact the case, not only outside the marginally 
stable orbit, but also inside, where heating by the returning radiation dominates. As a 
result, all the heat generated is released at the same radius, and the outcoming flux is 
the sum of the profiles of viscous and returning radiation heating rates (the blue solid 
line follows the standard solution outside the marginally stable orbit, and the returning 
radiation profile inside this radius). The maximal disk temperature, which determines 
the spectral hardness and is crucial for estimating the BH spin using X-ray continuum 
fitting (see Section 8.2), does not change — the self- irradiation affects the emission pro- 
file only inside the marginally stable orbit and the modified emission does not exceed 
the maximal value reached at r ^ lOM. 

The behavior discussed above changes when the accretion rate becomes high. For 
l.OMEdd the disk is no longer thin {h/r ^ 0.3), and the profile of emission of the 
standard solution slightly shifts inward due to advection. Furthermore, the disk surface 
is no longer straight (6* = const), but the innermost region is puffed up by the pressure 
of radiation. As a result, the inner region of this radiation pressure-related hump 
captures back radiation emitted from disk inner regions on the other side of the central 
BH. What follows, the radiation from the innermost disk does not reach the most 
distant (r > lOOM) regions, being eclipsed by the inner hump. The increase of the 
amount of returning radiation due to the radiation pressure-induced disk thickening 
around r = 30M is visible in Fig. 7.2. The orange dashed line presents its profile for 
l.OMEdd- At large distances it contributes to ~ 3% of the viscous heating, exactly as in 
the case of a thin disk. In the region inside r = 30M, however, its amount, compared to 
the standard solution, increases to ~ 20%, due to the photon capture by the hump. The 
profile of the returning fiux further increases towards the BH horizon due to efficient 
light bending close to the BH. All the captured radiation, as will be discussed in one of 
the following paragraphs is reemitted, and the profile of emission of the self-irradiated 
solution is, similarly as for the O.lMEdd case, the sum of the viscosity and returning 
fiuxes. This time, however, the maximal effective temperature of the emitted flux is 
affected and increases by ~ 5% when compared to the standard solution. 

The impact of the radiation pressure-induced inner hump, and, in general, increased 
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disk thickness is even more profound for accretion disks with super-critical accretion 
rates. The profile of the estimated location of the disk photosphere for lO.OMEdd is 
presented in Fig. 7.1. It is clear that disk inner region irradiates itself, and only small 
fraction of radiation is able to escape to the observer at infinity. The photons emitted 
from the most energetic, inner regions are not able to reach the most distant regions, 
obscured by the thickening. The green curves in Fig. 7.2 show this behavior in terms 
of the returning and emitted radiation. The flux re-captured by the disk in the inner 
region (r < 50M) is equal to ~ 80% of the emerging one — only fifth part of the emitted 
radiation escapes to infinity. If the disk was radiatively efficient, this flux of returning 
radiation would increase the emitted one almost by a factor of two. However, at such 
high accretion rate the disk is advection dominated, and most of the heat returning 
with radiation to the disk is advected towards BH instead of being instantaneously 
reemitted. As a result, the flux emitted by the self- irradiated solution is higher than 
the flux of the standard one only by ~ 25%, what corresponds to the same increase of 
the effective temperature as in the case of l.OMEdd- 

7.3.2 Returning radiation 

In Fig. 7.3 we plot the radial distribution of the radiation returning back to the disk at 
two radii: r = lOM (solid) and r = lOOM (dashed lines) for the same three accretion 
rates. The blue curves correspond to the lowest, O.lMEdd- The profiles of the returning 
fiux for both radii refiect the profile of the emitted radiation with the peak around lOM 
and a minimum slightly inside the marginally stable orbit. The returning fiux, however, 
does not increase, like the emitted fiux does, towards the horizon, as the photons 
coming back from the vicinity of BH are subject to very strong redshift making their 
contribution to the returning fiux at larger radii negligible. Even for large distances 
from the BH (e.g., r = lOOM, blue dashed line), the photons from the most energetic, 
inner regions around r = lOM contribute to the majority of the returning fiux. 

Once the disk thickens and the inner hump appears, the radial distribution of the 
returning radiation changes. Most of the fiux returning to r = lOM still comes from the 
hottest region inside this radius, the radiation from the vicinity of BH is also strongly 
redshifted and does not contribute significantly to the returning fiux. However, the 
distribution does no longer extend uniformly towards large radii. The region outside 
r 80M is obscured by the inner hump and no single photon from this region is 
able to return to the disk at r = lOM. This fact is clearly visible also on the profiles 
of the radiation returning to r = lOOM (dashed lines). For l.OMEdd the distribution 
terminates at r = lOM meaning that the region inside this radius is hidden for an 
observer located at the disk photosphere at r = lOOM. For lO.OMEdd this obscuration 
takes place starting from r = 20M. 

7.3.3 Energy balance 

At each radius the balance between heating and cooling must be satisfied. In the stan- 
dard approach, only three mechanisms were involved: viscous heating, radiative cooling 
and advective heating (decreasing radial advective fiux heats the disk at given radius 
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Figure 7.3: Radial profiles of radiation coming back to r = lOM (solid) and r = lOOM 
(dashed lines). The curves show where the flux returning to a given radius comes from. 
Their integrals are equal to the returning radiation flux at these radii plotted in Fig. 7.2. 
Profiles for three accretion rates: 0.1 (blue), 1.0 (orange) and lO.OMEdd (green lines) 
are presented. 
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up) or cooling (some fraction of heat generated by viscosity increases the advective 
flux). When the self-irradiation is taken into account, an additional heating mechanism 
resulting from the returning radiation must be considered. 

In Fig. 7.4 we compare the heating/cooling profiles of the standard (dashed) and 
self-irradiated (solid lines) slim disk solutions for three accretion rates: 0.1 (top), 1.0 
(middle) and lO.OMEdd (bottom panel). Advection for the lowest accretion rates (top 
panel), corresponding to radiatively efficient disks, is negligible. Therefore, in the stan- 
dard model the heat generated by viscosity is completely balanced by the radiative 
cooling. When the returning radiation is considered, another heating mechanism oc- 
curs. The excess of heating at r < 6M must be balanced by the radiative emission. 
This is in fact the case — the rate of heating by the returning radiation (brown solid 
line) is equal to the absolute value of the radiative cooling rate in this region (blue solid 
line) . 

Advection becomes an important factor in the energy balance when the accretion 
rate is large. The middle panel of Fig. 7.4 presents the heating and cooling rates 
for both models for l.OMEdd- In the standard approach, viscosity is the only heating 
mechanism at r > lOM. Small fraction (compare Fig. 4.10) of the heat generated in this 
region is stored in the advective flux and is released, contributing to the total heating 
rate, inside this radius. At 5M < r < 7M the advective heating dominates over the 
viscosity-induced mechanism, leading to slight inward shift of the emission. 

This characteristic changes, mostly in the innermost region, when the returning 
radiation is considered (solid lines) . The heating caused by this process is most efficient 
near the BH, where it dominates over the other mechanisms. For the whole region 
inside r = lOM, the advection contributes to heating only, and all the heat coming with 
the returning radiation must be balanced by the increase of radiative cooling. As a 
result, the amount of the emitted flux increases by the rate of the captured radiation, 
leading to the modified emission profile presented in Fig. 7.2. 

The solutions for the highest accretion rates (e.g., lO.OMEdd, bottom panel) are 
significantly different, as the advective cooling is prominent at most radii. For the 
standard model, the advection contributes to cooling at r > 4M and may be twice 
as efficient as the radiative mechanism at moderate (r ^ lOM) radii. This leads to 
an inefficient emission (blue line) equal to only a small fraction of the heating rate 
(orange). Disks at such high accretion rates are very thick and capture significant part 
of the emitted radiation. However, most of this radiation is not radiated away at the 
same radius due to the efficient advective cooling mechanism. Instead, the captured 
heat is advected towards the BH. As a result, the emission profile does not change 
significantly. Only the rate of advective cooling significantly increases. 

7.4 Discussion 

We have introduced a model of an advective, optically thick accretion disk which prop- 
erly accounts for the heating impact of the returning radiation, as well as the angular 
momentum taken away and brought in by photons. In the previous section we de- 
scribed in detail solutions for disks around non-rotating BHs and compared them with 
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Figure 7.4: Energy balance as a function of radius for three accretion rates: 0.1 (top), 
1.0 (middle) and lO.OMgdd (bottom panel). Positive values denote heating, while neg- 
ative stand for cooling. Profiles for the self-irradiated solutions are marked with solid 
lines, while dashed lines correspond to the regular disk model. Viscous heating, ra- 
diative cooling and advective transport are marked with orange, blue and green lines, 
respectively. The brown lines denote heating by the returning flux of radiation. 
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the solutions of the standard model presented in Chapter 4. This analysis leads to the 
following conclusions: 

(i) The impact of self-irradiation onto disks with sub-Eddington accretion rates is 
negligible, as the returning radiation is comparable with the viscous heating only inside 
the marginally stable orbit. Despite the fact that the emission profile is significantly 
modified in this region, the resulting spectrum will be hardly affected, as the radiation 
coming from the vicinity of BH is strongly redshifted and hardly contributes to the 
total emission. 

(ii) For disks with accretion rates close to the Eddington value, the relative amount 
of the captured radiation in the innermost region increases due to radiation pressure- 
induced disk thickening. For the Eddington accretion rate, the returning radiation 
around r = lOM may be equal up to 20% of the flux emitted there. At these accretion 
rates, all the returning radiation contributes to the emission, increasing the amount 
of the emerging flux. As a result, the maximal effective temperature increases by a 
few percent, leading to a spectrum slightly harder than the predicted by the standard 
model. 

(iii) The relative amount of the captured radiation is much larger for super-critical 
accretion rates. For lOMEdd, fhe returning flux may reach as much as 80% of the 
emitted one. However, at this regime of luminosities accretion disks are efficiently 
cooled by advection, and most of this returning heat is advected onto the BH. Only 
small fraction (20% for lOMEdd) contributes to the radiative emission. The emerging 
spectrum is slightly hardened in comparison with the standard model. 

The study presented in this Chapter was limited to the case of a non-rotating BH. 
In a forthcoming paper we will extend this approach to study the impact of returning 
radiation on accretion disks around rotating BHs. It has been shown by Ll ET AL. 
(2005) that its impact on disk spectra may be more significant for rapidly rotating 
BHs^. Our results show, that for slowly rotating central objects, only solutions with 
low accretion rates are unaffected by the returning radiation, and the standard models 
may be applied to model spectra. For higher accretion rates, however, they have to be 
used with CdbTQ^ clS they underestimate the amount of emitted radiation and the spectrum 
hardness. This discrepancy may be more profound for rotating BHs. Nevertheless, the 
uncertainty related to the model inadequacy is most likely negligible when compared 
to other unknown factors, e.g., the value of the viscosity parameter a at different 
luminosities (see Section 8.2). 



^They prove, that for the thin disk case, the impact of self-irradiation on the spectrum may be 
simulated to quite high accuracy by a slight increase of the accretion rate. 
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Chapter 8 
Applications 



8.1 Spinning- up black holes with super-critical accre- 
tion flows 

Astrophysical BHs are very simple objects - they can be described by just two parame- 
ters: mass and angular momentum. In isolation, BHs conserve the birth values of these 
parameters but often, e.g., in close binaries or in active nuclei, they are surrounded by 
accretion disks and their mass and angular momentum change. Accretion of matter 
always increases BH's irreducible mass and may change its angular momentum, usually 
described by the dimensionless spin parameter a* = a/M. The sign of this change and 
its value depend on the (relative) sign of accreted angular momentum and the balance 
between accretion of matter and various processes extracting BH's rotational energy 
and angular momentum. 

The question about the maximal possible spin of an object represented by the Kerr 
solution of the Einstein equation is of fundamental and practical (observational) inter- 
est. First, a spin a* > 1 corresponds to a naked singularity and not to a BH. According 
to the Penrose cosmic censorship conjecture, naked singularities cannot form through 
actual physical processes, i.e. singularities in the Universe (except for the initial one in 
the Big Bang) are always surrounded by event horizons (Wald, 1984). This hypothesis 
has yet to be proven. 

In any case, the "third law" of BH thermodynamics (Bardeen ET AL., 1973) 
asserts that a BH cannot be spun-up in a finite time to the extreme spin value a^^ = 1. 
Determining the maximum value of BH spin is also of practical interest because the 
efficiency of accretion through a disk depends on the BH's spin value. For example, 
for the "canonical" value a^. = 0.998 (see below) it is about fj ^ 32%, while for a* — 1 
one has fj -)■ 42% (Section 2.3). Recently, Banados ET AL. (2009) showed that the 
energy Eqom of the center-of-mass collision of two particles colliding arbitrary close to 
the BH horizon, grows to infinity {Eqom oo) when a* — 1^. 

A definitive study of the BH spin evolution will only be possible when reliable, non- 

^Of more fundamental interest is the fact that the proper geodesic distances D between the ISCO 
and several other special Keplerian orbits relevant to accretion disk structure tend to infinity D ^ oo 
when a* —J' 1 (Bardeen et al., 1972). 
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stationary models of accretion disks and jet emission mechanisms are established. For 
now, one has to use simplified analytical or numerical models. 

Thorne (1974) used the model of a radiatively efficient, geometrically thin ac- 
cretion disk (NoviKOV & Thorne, 1973) to evaluate BH spin evolution taking into 
account the decelerating impact of disk-emitted photons. The maximum value so ob- 
tained, a* = 0.9978, has been regarded as the canonical value for the maximal BH spin. 
In this Chapter we generalize Thome's approach, using models of advective, optically 
thick accretion disks (slim disks. Chapter 4) to calculate maximum BH spin values for 
a large range of accretion rates. We show that for sufficiently large accretion rates they 
differ from the canonical value. 

We start with a short discussion of previous works devoted to the BH spin evolution. 
In Section 8.1.2 we give formulae for a general tetrad of an observer comoving with the 
accreting gas along the arbitrary photosphere surface. In Section 8.1.3 we give basic 
equations describing the BH spin evolution. In Section 8.1.4 we present and discuss 
the terminal spin values for all the models we consider. Finally, in Section 8.1.5 we 
summarize our results. 

8.1.1 Previous studies 

A number of authors have studied the evolution of the BH spin resulting from disk 
accretion. Bardeen (1970) initiated this field of research by stating the problem and 
solving equations describing the BH spin evolution for accretion from the marginally 
stable orbit. Neglecting the effects of radiation he proved that such a process may 
spin-up the BH up to a* = 1. Once the classical models of accretion disks were for- 
mulated (Shakura & SuNYAEV, 1973; NoviKOV & Thorne, 1973), it was possible 
to account properly for the decelerating impact of radiation (frame dragging makes 
counter-rotating photons more likely to be captured by the BH). As mentioned above, 
Thorne (1974) performed this study and obtained the terminal spin value for an 
isotropically emitting thin disk of a* = 0.9978, independently of the accretion rate. 
The original study by Thorne was followed by many papers, some of which are briefly 
mentioned below. 

The flrst to challenge the universality of Thome's limit were AbramowiCZ & LA- 
SOTA (1980) who showed that geometrically thick accretion disks may spin up BHs to 
terminal spin values much closer to unity than the presumably canonical a^, = 0.9978. 
Their simple argument was based on models by KOZLOWSKI ET AL. (1978) who showed 
that for high accretion rates the inner edge of a disk may be located inside the marginally 
stable orbit; with increasing accretion rate, arbitrarily close to the marginally bound 
orbit. However, this conclusion assumed implicitly a low viscosity parameter a, whereas 
for high viscosities the situation is more complicated (see Abramowicz ET AL. (2010) 
and references therein). 

MODERSKI ET AL. (1998) assessed the impact of possible interaction between the 
disk magnetic fleld and the BH through the B landlord- Znajek process. They showed 
that the terminal spin value may be decreased to any, arbitrarily small value, if only 
the disk magnetic fleld is strong enough. Given the current lack of knowledge about 
the strength of magnetic flelds (or the magnetic transport of angular momentum in the 
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disk, see Ghosh & Abramowicz (1997)) and processes leading to jet emission, a 
more detailed study cannot be performed. The situation may further be complicated 
by energy extraction from the inner parts of accretion disks (LiVIO ET AL., 1999). 

POPHAM & Gammie (1998) studied the spinning-up of BHs by optically thin advec- 
tion dominated accretion flows (ADAFs). They neglected the contribution of radiation 
to BH spin as such disks are radiatively inefficient. They found that the terminal value 
of BH spin is very sensitive to the assumed value of the viscosity parameter a and may 
vary between 0.8 and 1.0. Gammie et AL. (2004), besides making a comprehensive 
summary of different ways of spinning up supermassive BHs, presented results based 
on a single run of a general relativistic magnetohydrodynamical (GRMHD) simulation 
(with no radiation included) obtaining terminal spin a* = 0.93. 

Cosmological evolution of spins of supermassive BHs due to hierarchical mergers and 
thin-disk accretion episodes has been recently intensively studied. Although VOLON- 
TERI ET AL. (2005) arrived to the conclusion that accretion tend to spin-up BHs close 
to unity, as opposed to mergers which, on the average, do not influence the spin evolu- 
tion, the following studies by e.g., VOLONTERI ET AL. (2007); KiNG ET AL. (2008); 
Berti & VOLONTERI (2008) showed that the situation is more complex, the final 
spin values depending on the details of the history of the accretion events (see also 
Fanidakis et AL. (2011)). 

Li et AL. (2005) included the returning radiation into the thin-disk model of 
NoviKOV & Thorne (1973) and calculated the spin-up limit for the BH assuming the 
radiation crossing the equatorial plane inside the marginally stable orbit to be advected 
onto the BH. Their result (a* = 0.9983) slightly differs from Thome's result, thus 
showing that returning radiation has only a slight impact on the process of spinning-up 
BHs. In our study we use advective, optically thick solutions of accretion disks and 
account for photons captured by the BH in detail. However, we neglect the impact of 
the radiation returning to the disk on its structure. 

8.1.2 Tetrad 

We base this work on slim accretion disks, which are not razor-thin and have angular 
momentum profile that is not Keplerian (for details on the assumptions made and the 
disk appearance see Chapter 4). Therefore, photons are not emitted from matter in 
Keplerian orbits in the equatorial plane and the classical expressions for photon mo- 
menta (e.g., MiSNER ET AL. (1973)) cannot be applied. Instead, to properly describe 
the momentum components of emitted photons, we need a tetrad for the comoving 
observer instantaneously located at the disk photosphere. Below we give the explicit 
expression for the components of such a tetrad assuming time and axis symmetries. A 
detailed derivation is given in Appendix C. 
Let us choose the following comoving tetrad, 

e;^) = K,iV:,«:^,5^ (8.1) 

where is the four- velocity of matter (living in [t, </>, r, 6*]), 
A^* is a unit vector orthogonal to the photosphere ( [r, Q\ ) , 
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unit vector orthogonal to ([t, </>]), 

is a unit vector orthogonal to m*, A^* and Kq ([t, </>, r, 6]). 
The tetrad components are given by, 



-1/2 



d^ 
dr 



5'r 



V dr 



1 -1/2 



-1/2 



+ fie + ^''5: 



1.5) 



S' = {1 + A^vY'^^iAvu' + Si), 



].6) 



where 6 = 9^{r) defines the location of the photosphere, rji and ^ are the Killing vectors, 
/ = u^/ut, fi = u'^/u'^, u is the frequency of frame-dragging and the expressions for SI 
and V are given in Eqs. C.5 and C.ll, respectively. 



8.1.3 Spin evolution 

Basic equations 

The equations describing the evolution of BH dimensionless spin parameter a* with 
respect to the BH energy M and the accreted rest-mass Mq are (Thorne, 1974), 
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rad 



dlnM dlnM M MoUt + {^) 



2a*, 



rad 



dM fdM\ 

Ut+{ — ] /Mo. 



dMo 



dt 



rad 



The energy and angular momentum of BH increases due to the capture of photons 
according to the following formulae. 
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where rjk and are the Kilhng vectors connected with time and axial symmetries, re- 
spectively, T*'^ is the stress-energy tensor of photons, which is non-zero only for photons 
crossing the BH horizon and dS, the "volume element" in the hypersurface orthogonal 
to A^*, is given by Eq. D.8. 

From Eqs. 8.9 and 8.10 it follows that 

i^j = I I r%NdS, (8.11) 



rad JO 



2tT /Tout 



where 



dS = d0dr [gl^ - gug^p^Y^"^ J gw + gee {^^^ ■ (8-13) 



Stress energy tensor in the comoving frame 

Let us choose the tetrad given in Eq. 8.1: 

e%^=u^ e\,) = N^ (8.14) 
6(2) = ejg) = S' 

The disk properties, e.g., the emitted flux, are usually given in the comoving frame 
defined by Eq. 8.14. The stress tensor components in the two frames (Boyer-Lindquist 
and comoving) are related in the following way, 

= T(")(^)eJ„)ef^). (8.15) 
The stress tensor in the comoving frame is 

/'7r/2 p2n 

rp{am = 2 / loSCn^^h^'^hmadddb, (8.16) 

Jo Jo 

where I^S = Io{r)S{a, b) is the intensity of the emitted radiation, a and b are the angles 
between the emission vector and the N'^ and S** vectors, respectively, C is the capture 
function defined in Section 8.1.3, the factor 2 comes from the fact that the disk emission 
comes from both sides of the disk and are the normalized components 

of the photon four-momentum in the comoving frame. The latter are given by the 
following simple relations (Thorne, 1974), 

7r(°) = 1, 

7r« = cos a, (8.17) 

TT*-^^ = sin a COS 6, 

TT*-^-' = sin a sin 6. 
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Eqs. 8.11 and 8.12 take the form 

'dM 



dt 



T(")('^)e;„)ef^)77,iV,,d^, (8.18) 

rad ^disk 



(^) = [ T(")(^)e;,)ef,)aiV.d^, (8.19) 

V / rad ^disk 

where e^^-j is our local frame tetrad given by Eq. (8.14). Taking the following relations 
into account, 

= ^igW^ (820) 



we have, 



where, 



7r(")e;,)iV, = 7r(")5j;;5 = tt^^) = cos a, (8.21) 

7r(^)ef^)% = Ti'efu\p^7]u = n'S^r], = n% = Tit, (8.22) 

vr^^^ef^)^ = vr^ef ^wf^)^ = vr^^J^ = ^'^k = vr^, (8.23) 

TTf = 7r«e*,)^?« + 7r«eJ)^?i^, (8.24) 

TT^ = 7r«eJ,)(7t<^ + 7r«eJ)(7^<^. (8.25) 

Therefore, Eqs. 8.18 and 8.19 may be finally expressed as, 

' ^ = in / / loSCnt X 

rad ^r-in ^0 JO 



V dt 



X 



cos a sm a dadb^dr, (8.26) 



J T\ rrout r2n /•7r/2 

= 47r / / / IoSCtt^ X 

rad -'nn >^0 Jo 



dt 



X 



COS a sm a dddb^/^dr, (8.27) 



with 



^ 9tt9<p^Y^'^ \l 9rr + gee (^j ■ (8.28) 
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Emission 

The intensity of local radiation may be identified with the flux emerging from the disk 
surface F{r) (presented and discussed for the case of slim disks in Section 4.3) 

/o = F(r). (8.29) 

The angular emission factor S is given by (Thorne, 1974), 

S(a b) = I ^^'^ isotropic i'^'i^i) 

(3/77r)(l + 2 cos a) limb darkening 

for isotropic and limb-darkened cases, respectively. In this work we assume that the 
radiation is emitted isotropically. 

Capture function 

The BH energy and angular momentum are affected only by photons crossing the BH 
horizon. Following Thorne (1974), we define the capture function C, 

^ _ r 1 if the photon hits the BH, ,^ 

in the opposite case. ^ ' ' 

Herein, we calculate C in two ways. First, we use the original Thorne (1974) algo- 
rithm modified to account for emission out of the equatorial plane. For this purpose 
we calculate the constants of motion, j and /c, for a geodesic orbit of a photon in the 
following way, 

J = a' + a*(WM7rt), (8.32) 

fc = VI - + sin Q,fl sin^ Q,\ (8.33) 

which replaces Thome's Eqs. AlO. This approach does not take into account possible 
returning radiation, i.e., a photon hitting the disk surface is assumed to continue its 
motion. Such a treatment is not appropriate for optically thick disks - returning photons 
are most likely absorbed or advected towards the BH. 

To assess the importance of this inconsistency we adopt two additional algorithms 
for calculating C . Using photon equations of motion we determine if the photon hits 
the disk surface (BURSA, 2006). Then we make two assumptions, either the angular 
momentum and energy of all returning photons are advected onto the BH (Ci), or all 
are reemitted carrying away their original angular momentum and energy, and never 
hit the BH (C2). In this way we establish two limiting cases allowing us to assess the 
impact of the returning radiation. 

We note here that for fully consistent treatment of the returning radiation (as in 
Ll ET AL. (2005) for geometrically thin disks) it is not enough to modify the capture 
function — solving for the whole structure of a self-irradiated accretion disk (like in 
Chapter 7) is necessary. The impact of self-irradiation onto the BH spin evolution will 
be studied in detail in the future. 
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8.1.4 Results for BH spin evolution 

Using the slim disk solutions described in Chapter 4 we solve Eqs. 8.7 and 8.8 using 
regular Runge-Kutta method of the 4th order. To calculate the integrals (Eqs. 8.11 
and 8.12) we use the alternative extended Simpson's rule (PRESS, 2002) basing on 100 
grid points in a, b and radius r. We have made convergence tests proving this number 
is sufficient. 

In Figs. 8.1 and 8.2 we present the BH spin evolution for a = 0.01 and 0.1, respec- 
tively. The red lines show the results for different accretion rates while the black line 
shows the classical Thorne (1974) solution based on the NoviKOV & Thorne (1973) 
model of thin accretion disk. The accretion rates are given in the units of the critical 
accretion rate defined in Eq. 4.35. Our low accretion rate limit does not perfectly agree 
with the black line as the slim disk model does not account for the angular momentum 
carried away by radiation. As a result, the low-luminosity slim disk solutions slightly 
overestimate (by no more than few percent) the emitted flux leading to stronger decel- 
eration of the BH by radiation — the Thorne (1974) result is the proper limit for 
the lowest accretion rates. When the accretion rate is high enough (e.g., m > 0.1), the 
impact of the omitted angular momentum flux is overwhelmed by the modification of 
the disk structure introduced by advection. 

It is clear the spin evolution is quantitatively different for different values of the 
viscosity parameter. For the lower value {a = 0.01) the BH spin can reach values 
significantly higher than 0.998 for the highest accretion rates {a^ = 0.9995 for fn = 100), 
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Figure 8.2: Spin evolution for a = 0.1. 



while for the higher viscosity (a = 0.1) the terminal spin value decreases with increasing 
accretion rate down to = 0.9800 for rh = 100. This behavior is connected with the 
impact of viscosity on the angular momentum profile. Generally speaking, the higher 
the value of a, the lower the angular momentum of the flow at BH horizon for a given 
accretion rate (see Fig. 4.8 in Section 4.3), leading to a slower acceleration of the BH 
rotation. 

To study this fact in detail we calculated the rate of BH spin-up for "pure" accretion 
of matter (without accounting for the impact of radiation). For that case the BH spin 
evolution is given by (compare Eq. 8.7): 



- ^ - 2a*. (8.34) 



dlnM M ut 

In Figs. 8.3 and 8.4 we plot with black lines the first term on the right hand side 
of the above equation for different accretion rates and values of a. The red lines on 
these plots show the absolute value of the second term. The intersections of the black 
and red lines denote the equilibrium states i.e., the limiting values of BH spin for pure 
accretion. These values differ significantly from the previously discussed results only for 
low accretion rates. In contrast, for high accretion rates radiation has little impact on 
the spin evolution and the value of terminal spin is mostly determined by the properties 
of the flow. In Fig. 8.5 we plot the radiation impact parameter ^, defined as the ratio 
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Figure 8.3: The rate of spin-up or spin-down by "pure" accretion (radiation neglected) 
for a = 0.01. Profiles for five accretion rates are presented. Their intersections with 
the red line (marked with blue crosses) correspond to equilibrium states. For the two 
lowest accretion rates the equilibrium state is never reached (a* — )■ 1). 

of the disk-driven terms on the right hand sides of Eqs. 8.7 and 8.34, 

If the captured radiation significantly decelerates BH spin-up, this ratio drops below 
unity. On the other hand, it is close to unity for BH spin evolution which is not 
affected by the radiation. According to Fig. 8.5 the latter is in fact the case for the 
highest accretion rates independently of a. 

In Table 8.1 we list the resulting values of the terminal BH spin for all the models 
considered. The first column gives results for our fiducial model (A) including Thome's 
capture function, emission from the photosphere at the appropriate radial velocity. 

The second column presents results obtained assuming the same (Thome's) capture 
function and profiles of emission, angular momentum and radial velocity as in the 
model A, but assuming the emission takes place from the equatorial plane instead of 
the photosphere. The resulting terminal spin values are equal, up to 4 decimal digits, 
to the values obtained with the fiducial model. This result is expected for the lowest 
accretion rates, where the photosphere is located very close to the equatorial plane. 
For the highest accretion rates the location of the emission has no impact on the BH 
spin-up, as the spin evolution is driven by the flow itself and the effects of radiation are 
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Figure 8.4: Same as Fig. 8.3 but for a = 0.1. 



negligible. However, for moderate accretion rates one could have expected significant 
change of the terminal spin. We find that the location of the photosphere has little 
impact on the resulting BH spin regardless of the accretion rate. 

Our third model (V) neglects the flow radial velocity when the radiative terms are 
evaluated. Similar arguments to those given in the previous paragraph apply. For 
the lowest accretion rates, the radial velocity is negligible and therefore it should have 
no impact on the resulting spin. For the highest accretion rates, the spin-up process 
depends only on the properties of the flow. Once again, however, the impact of this 
assumption on moderate accretion rates is not obvious. The radial velocity turns out 
to be of little importance for the calculation of the terminal spin (only for a = 0.1 and 
moderate accretion rates the difference between models A and V is higher than 0.01%). 

In the fourth and fifth columns of Table 8.1, results for models with the same as- 
sumptions as the fiducial model, but with different capture functions are presented. The 
first alternative capture function (Ci) assumes that the angular momentum and energy 
of all photons returning to the disk are added to that of the BH. This assumption has a 
strong impact on the spin evolution — the terminal spin values are higher, sometimes 
approaching a* = 1. This may seem surprising because in the classical approach the 
captured photons are responsible for decelerating the spin-up. This effect comes from 
the fact that the cross-section (with respect to the BH) of photons going "against" the 
frame dragging is larger than of photons following the BH sense of rotation. As frame 
dragging is involved, this effect is significant only in the vicinity of the BH horizon. 
For our model Ci, however, the probability of photons returning to the disk does not 
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Figure 8.5: Radiation impact factor ^ (Eq. 8.35) for different accretion rates and values 
of a. The dotted line corresponds to the thin-disk induced spin evolution. For ^ ^ 1 
spin evolution is not affected by radiation. Stars denote the equilibrium states (compare 
Table 8.1). 
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Table 8.1: BH spin terminal values. C - Thome's capture function, Ci - all returning 
photons advected onto the BH, C2 - all returning photons neglected; A - our fiducial 
model, T - emission from the equatorial plane, V - zero radial velocity, NR - pure 
accretion, radiation neglected. 

appreciably differ for co- and counter-rotating photons, as they both hit the disk surface 
mostly at large radii. 

The other capture function (C2) assumes, on the contrary, that all returning photons 
are re-emitted from the disk with their original angular momentum and energy (and 
never fall onto the BH). This assumption cuts off the photons which would hit the BH 
in the fiducial model after crossing the disk surface. Thus, we may expect decreased 
radiative deceleration and increased values of the terminal spin parameter. However, 
these changes are not significant, refiecting the fact that most of the photons hit the BH 
directly, along slightly curved trajectories. Only for a = 0.1 and moderate accretion 
rates do the terminal spin values differ in the 4**^ decimal digit. 

Neither of the models with a modified capture functions is self-consistent. To ac- 
count properly for the returning radiation one has to modify the disk equations by 
introducing appropriate terms for the outgoing and incoming fiuxes of angular momen- 
tum and additional radiative heating. No such model for advective, optically thick 
accretion disk has been constructed. The emission profile should be significantly af- 
fected (especially inside the marginally stable orbit) by the returning radiation leading 
to different rates of deceleration by photons. In view of our results for models Ci and 
C2, as well as the results of Ll ET AL. (2005), one may expect the final spin values for 
super-critical accretion fiows to be slightly higher than the ones obtained in this work. 

The last column of Table 8.1 gives terminal spin values for "pure" accretion (radiation 
neglected) . Under such assumptions the BH spin could reach a* = 1 for sub-Eddington 
accretion rates as there are no photons which could decelerate and stop the spin-up 
process. As discussed above, for the highest accretion rates the resulting BH spin values 
agree with the values obtained for the fiducial model as radiation has little impact on 
spin evolution in this regime. 
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8.1.5 Summary 

We have studied BH spin evolution due to disk accretion assuming that the angular 
momentum and energy carried by the flow and the emitted photons is the only process 
affecting the BH rotation. We generalized the original study of Thorne (1974) to 
high accretion rates by applying a relativistic, advective, optically thick slim accretion 
disk model. Assuming isotropic emission (no limb darkening) we have shown that 

(i) the terminal value of BH spin depends on the accretion rate for m > 1, 

(ii) the terminal spin value is very sensitive to the assumed value of the viscosity 
parameter a — for a < 0.01 the BH is spun up to a* > 0.9978 for high accretion rates, 
while for a > 0.1 to a, < 0.9978, 

(iii) with a low value of a and high accretion rates, the BH may be spun up to 
spins significantly higher than the canonical value a* = 0.9978 (e.g., to a* = 0.9994 
for a = 0.01 and rh = 10) but, under reasonable assumptions, BH cannot be spun up 
arbitrarily close to a^, = 1, 

(iv) BH spin evolution is hardly affected by the emitted radiation for high (m > 10) 
accretion rates (the terminal spin value is determined by the flow properties only), 

(v) for all accretion rates, neither the photosphere profile nor the profile of radial 
velocity significantly affects the spin evolution. 

We point out that the inner edge of an accretion disk cannot be uniquely defined 
for a super-critical accretion (Section 4.4), as opposed to geometrically thin disks where 
the inner edge is uniquely located at the marginally stable orbit (rigco)- In the thin 
disk case the BH spin evolution is determined by the flow properties at this particular 
radius (as there is no torque between the marginally stable orbit and BH horizon) 
and the profile of emission (terminating at rigco)- For super-critical accretion rates, 
however, one cannot distinguish a particular inner edge which is relevant to studying 
BH spin evolution. On the one hand, the values of the specific energy (ut) and the 
angular momentum (u^) remain constant within the stress inner edge. On the other, 
the radiation is emitted outside the radiation inner edge. These inner edges do not 
coincide as they are related to different physical processes. 

Our study was based on a semi-analytical, hydrodynamical model of an accretion 
disk which makes a number of simplifying assumptions like stationarity, no returning 
radiation, ap prescription and no wind outflows. Thus, the results obtained in this 
work should be considered only qualitative, as the precise values of the terminal spin 
parameter are very sensitive (e.g., through viscosity) to the flow and emission properties. 
Moreover, we neglected the impact of magnetic fields and jet ejection mechanisms. 
Nevertheless, our study shows that Thome's canonical value for BH spin {a^ = 0.9978) 
may be exceeded under certain conditions. 

8.2 Spin determination via X-ray continuum fitting 

Astrophysical BHs may be uniquely characterized by two parameters: mass and angular 
momentum. From among those two, the latter is much more difficult to measure, as it 
has no impact onto the orbital motion of the compact binary. Currently, three methods 
of estimating BH spin are in use: based on the shape of emission lines (e.g., Reynolds 
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& Fabian (2008)), on the continuum fitting (e.g., Steiner ET AL. (2010)), and on 
studying the resonances leading to quasi-periodic oscillations (QPOs, TOROK ET AL. 
(2007)). First two give consistent results for those few objects that can be studied 
using both approaches. 

The continuum-fitting method is usually based on the NOVIKOV & Thorne (1973) 
thin disk model and is used for X-ray binaries which exhibit clear high/soft spectra 
dominated by the disk thermal component (Section 1.4). It has been successfully applied 
for a number of objects. However, its range of applicability seems to be limited to low- 
luminosity spectra only. As has been shown (e.g., Steiner ET AL. (2010)), the spin 
estimates for high luminosity spectra are not consistent. The characteristic tail of 
dropping spin estimates in this regime (Fig. 8.6) is common for all objects. 

This inconsistency emerges for luminosities which are close to or outside the limit of 
applicability of radiatively efficient disks. We have shown that advective significantly 
modifies the disk structure for L > 0.6LEdd (Section 4.4). Therefore, thin disk models 
should not be applied in this regime, as they cannot produce consistent disk fiux profiles 
which are used for continuum fitting purposes. 

In this section we briefiy discuss the application of the spectral model slimbb (Sec- 
tion 4.5.2) to high/soft state observations of LMC X-3. We report results described in 
StRAUB ET AL. (2011). 

The data which has been analyzed, consisting of 712 PCU-2 X-ray spectra, was 
obtained using the large-area Proportional Counter Array (PCA) aboard the Rossi 
X-ray Timing Explorer (RXTE). The most important selection criteria for those set 
of spectra were: (i) disk component exceeding 75% of the total emission, (ii) small 
variability, (iii) no high-frequency QPOs^. The X-ray spectral fitting software, XSPEC 
(Arnaud, 1996) was used with the following model combination TBabs * simpl * 
slimbb, where TBabs stands for the Galactic neutral hydrogen absorption with uh = 
4 X 10^° cm"^ (Page ET AL., 2003) and simpl for a power-law component (Steiner 
ET AL., 2009). The mass M = IOMq, the average distance to the LMC, D = 48.1 kpc 
(OrOSZ ET AL., 2009) and the inchnation, i = 66° (KuiPER ET AL., 1988) were 
assumed. The proper shape of the Compton-hardened spectrum was taken from the 
BHSPEC atmospheric models (Davis et AL., 2005). 

In Fig. 8.6 we present a comparison of the fittings obtained with slimbb with analo- 
gous study conducted with kerrbb2 — the BHSPEC enhanced version of kerrbb (with 
switched off returning radiation, to match the slim disk assumptions). The top panel 
presents results for a = 0.1. Both approaches give the same results producing the 
characteristic, artificial spin-luminosity dependence — the results are in fact indistin- 
guishable. Only for low luminosities (L < 0.2LEdd) one gets roughly constant value of 
BH spin a* ~ 0.6. For higher luminosities the spin estimates drop down — as far as to 
a* = 0.3 for 0.5LEdd- 

The bottom panel presents the corresponding comparison of slimbb and kerrbb2 
models for lower value of the viscosity parameter a = 0.01. The fittings behave quali- 
tatively in the same way. However, the spin-luminosity relation is much fiatter — the 



■^QPOs are observed only in the low/hard state or during state transitions, so the data containing 
QPOs is most hkely not disk-dominated. 
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spin estimates drop down only to a* = 0.4 for O.SLEdd- Similarly as in the a = 0.1 case, 
the results provided by these models are indistinguishable. 

In comparison to the standard (based on thin disk models) kerrbb2, the model 
based on slim disk solutions (slimbb) accounts additionally for the advective cooling. 
In this way it is able to consistently describe accretion flows with moderate and high 
luminosities. As the comparison of BH spin estimates for LMC X-3 shows, this single 
improvement is not able to solve the high-luminosity inconsistency. Despite the fact 
that advection modifies the emission (however, as discussed in Section 4.3, its impact is 
hardly visible for sub-Eddington luminosities), the spectra produced by combined disk 
and atmospheric models are too hard, when compared to the observed data (the harder 
the spectrum, the lower the spin estimate). One may conclude that the reason for this 
discrepancy lies in a different process than advection. 

There are at least few shortcomings of the hydrodynamical models in use, which 
may be responsible for this fact. To start with, these models neglect mass outflows 
assuming constant (in radius) accretion rate. However, it has been recently shown that 
radiation driven outflows, while hardly modify the effective temperature, lead to signifi- 
cant decrease of the local surface density which, in turn, would lead to further spectrum 
hardening (TakeuCHI ET AL., 2009). Another possible source of this discrepancy may 
be connected to the treatment of comptonization in the disk atmosphere. Currently, 
the disk atmosphere model BHSPEC, based on radiative transfer code TLUSTY, is the 
only available. Its approach to comptonization is sophisticated, but still involves some 
approximations. It would be beneficial, if spectra produced by BHSPEC may be com- 
pared to another similar, stand-alone radiative transfer code. It is, most probably, the 
question of time. 

The viscosity in accretion disks most likely results from magnetorotational insta- 
bility (Section 2.2.2). Magnetic fields are not explicitly taken into account in hydro- 
dynamical models (e.g., thin or slim disks). One has to include the viscosity in an 
artificial way. It is usually done by assuming that the stress is proportional to pressure 
(Section 2.2.1). Recent MHD shearing-box studies tend to support this assumption for 
mean accretion flows. These local simulations, however, say little about the dependence 
of the a parameter on radius or luminosity. Only recently, shearing boxes with different 
values of radiation to gas pressure ratio have been studied (HiROSE ET AL., 2009b) 
proving that the value of a, defined as the ratio of stress to local pressure, may slightly 
depend on the luminosity (Fig. 8.7). 

One may infer about the a dependence on luminosity not only from numerical sim- 
ulations of MRI, but also from observations of accretion disks at high luminosities. As 
Fig. 8.6 shows, the spin estimates are a-degenerate at low luminosities (L < 0.2LEdd)- 
At higher accretion rates, however, they strongly depend on a — the higher the value 
of a, the lower the spin value. This behavior results from the fact that lower a corre- 
sponds to higher surface density at given radius, which, in turn, decreases the efficiency 
of spectrum hardening through comptonization. If one sets the BH spin to a constant 
value corresponding to the results of the low-luminosity data analysis, then fitting for 
the value of a which produces a consistent spin estimate, is possible. 

In Fig. 8.8 we present results of such a study for LMC X-3 assuming BH spin 
a* = 0.61 (the mean value for fittings based on the L < 0.2LEdd data). The large error 
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bars for low-luminosity data reflect the fact that spin estimates are a-degenerate in this 
regime. For L > O.SLEdd? however, the magnitude of viscosity may be determined. This 
study shows that high (a > 0.05) values of a are forbidden when the disk luminosity 
is not small. Therefore, either the effective a in LMC X-3 is independent of luminosity 
and small {a < 0.01), or its value is larger for low-luminosity states and decreases 
with increasing luminosity down to the allowed values in the high luminosity regime. 
It is worth noting, that MHD shearing-box simulations of radiation dominated disks 
(Fig. 8.7) suggest weak decrease of a with luminosity (from 0.03 down to 0.01) what is 
both in qualitative and quantitative agreement with the results previously discussed. 



8.3 Normalization of global MHD disk simulations 

The general relativistic MHD models that are currently in use (e.g., Penna ET AL. 
(2010)) often make use of dimensionless quantities and do not include detailed radia- 
tion transfer. Hence, there is no direct way of estimating the physical mass accretion 
rate (g/s) or the true radiative luminosity (erg/s) of these models. To estimate these 
quantities one may use an indirect method, in which the vertical thicknesses of the sim- 
ulated disks is compared against physical disk models that include radiation transfer 
and radiation pressure and solve for the vertical disk structure, e.g., two-dimensional 
slim disks introduced in Chapter 6. 

Our solutions of the two-dimensional slim disk model (Chapter 6) were recently 
used, together with a non-LTE, radiatively efficient model of DAVIS ET AL. (2005), to 
estimate the luminosities that the set of global MHD simulations presented in KULKA- 
RNI ET AL. (2010) corresponds to. In this section we briefly compare their profiles of 
the disk thickness with the predictions of our two-dimensional slim disk model. 

In Fig. 8.9 we plot the density-weighted disk thickness, defined as, 

, , r„°° pz dz 

\h\ = ^° . 8.36 
Jo /'dz 

as predicted by the slim disk model (red solid lines, corresponding to L = 0.3^0.8LEdd, 
bottom to top) for a non-rotating BH. These curves should be compared with the 
disk thickness in the GRMHD simulations (deep blue lines). Results for two different 
configurations of the initial magnetic filed, purely-toroidal and loop, are shown with 
solid and dotted lines, respectively. The slim disk photosphere profiles are also presented 
(light blue, broken lines). 

At radii close to the ISCO, the slim disk and BHSPEC models indicate that the 
disk thickness plunges to small values whereas the GRMHD simulation shows a much 
smaller decrease. We believe there are at least three reasons for this discrepancy: (i) 
the GRMHD simulations cool the gas by forcing it towards a constant entropy, which 
may not be justified in the plunging region; (ii) the simulated GRMHD disk includes 
magnetic fields whose pressure provides additional support in the vertical direction 
whereas the other two models do not; and (iii) the simulated disc begins to deviate 
from hydrostatic equilibrium as the radial velocity becomes large near the ISCO and 
the gas has less time to reach equilibrium, whereas the other models hardwire the 
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condition of hydrostatic equilibrium at all radii. We estimate that the last two are only 
important well inside the ISCO (compare Fig. 6.13). These are interesting issues which 
we hope to explore in the future. 

For the purpose of comparing the disk thickness profiles, the region of the simulation 
near the ISCO is simply ignored. The comparison is based on radii r = 9 ^ 15M, which 
are well outside the ISCO and determines the luminosities at which the slim disk model 
gives the same disk thickness as obtained in the simulated GRMHD disk. We see 
from Fig. 8.9 that the thickness measure \h\/r ~ 0.05 in the simulated GRMHD disk 
corresponds to L/L-^dd ~ 0.45 according to the slim disk model. Similar analysis for 
a = 0.7, 0.9 andO.98 shows that L/LEdd ~ 0.35, 0.5 andO.5 respectively. 

Performing such comparison KULKARNI ET AL. (2010) were able to draw the fol- 
lowing conclusion: the observational errors in current measurements of BH spin by the 
continuum-fitting method dominate over the errors incurred by using the idealized NT 
model for L/Lsdd < 0.3. Comparison of the results of GRMHD simulations of accretion 
disks with the slim disk model was a crucial point in justifying this result. 
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Figure 8.6: Results of BH spin estimation for LMC X-3 using the continuum fitting 
method. Disk models slimbb (red diamonds) and kerrbb2 (orange triangles) for vis- 
cosity parameters a = 0.1 (top) and a = 0.01 (bottom) are compared. In all cases we 
assume LMC X-3 contains a lOM© black hole, is located at a distance of 48.1kpc and 
seen at an inclination of 66" (figure after: Straub et AL. (2011)). 
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Figure 8.7: Measured values of the stress parameter a as a function of the time-averaged 
ratio of the box-averaged radiation pressure to the box-averaged gas pressure. The black 
points define a as the time-averaged ratio of the vertically averaged stress to the box- 
averaged total thermal pressure. The blue and green points define a in the same way 
except with total thermal pressure replaced by gas pressure and the geometric mean of 
gas and radiation pressure, respectively. Both horizontal and vertical error bars indicate 
one standard deviation in the time averages. The horizontal black line indicates the 
weighted mean of a for the total pressure stress prescription. Assuming that stress 
really scales with total thermal pressure, the green and blue curves show how a would 
then behave under the other stress prescriptions, structure parameters. Figure and 
caption after: HiROSE ET AL. (2009b) 
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Figure 8.8: Behavior of the a- viscosity parameter with luminosity for LMC X-3. Figure 
after: Straub ET AL. (2011). 



H/R comparison, a=0.0 



0.1 



1 1 1 

: 

• / y 


1 1 1 



: 



;>> ; 

^ , ' 


1 1 

photosphere 

— — — — 






■ I / '/' / / / 

■ ,f i/' / / / 

' \g 1,' / / / y 

1 •'«'} // / / 1 


1 1 1 


L/Ledd=0.3 ' 

L/Ledd=0.4 

L/Ledd=0.5 

L/Ledd=0.6 

L/Ledd=0.7 

L/Ledd=0.8 

GRMHD 

GRMHDloop 

1 1 



4 6 8 10 12 14 16 18 20 

R[M] 

Figure 8.9: Comparison of the density- weighted disk thickness predicted by the two- 
dimensional slim disk model (red lines) with the results of two GRMHD simulations for 
a non-rotating BH. Photosphere profiles of the slim disk solutions are presented with 
violet dotted lines. 
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Chapter 9 
Conclusions 



Accretion of matter onto compact objects is one of the most efficient ways of extracting 
gravitational energy. Most luminous objects in the Universe: Active Galactic Nuclei, 
quasars, X-ray binaries are fed by accretion. Physics of accretion has, therefore, been 
one of the primary interests of astrophysicists for many years. 

Currently, two approaches to modeling accretion disks are in use: (semi-)analytical 
and numerical. The former was pioneered by the seminal work of Shakura & SUN- 
YAEV (1973). In this approach, the most general equations of magnetohydrodynamics 
are reduced, thanks to a number of simplifying assumptions, to a set of ordinary dif- 
ferential or algebraic equations. The most obvious advantage of such treatment is the 
simplicity of the final equations and, what follows, small computational complexity. 
The slim disk models described in this thesis are based on such an approach. 

The other way of modeling accretion disks is based on solving the equations de- 
scribing them in the full, unreduced form, using sophisticated numerical methods. In 
this way one is able to follow the magnetorotational instability in detail, as well as to 
discard all of the simplifying, but only approximate, assumptions. Unfortunately, such 
an approach is extremely computational time-consuming and therefore, its applicability 
is limited. A global (3D), relativistic MHD simulation of a thin disk with simplified 
treatment of cooling, similar to the one described in Section 8.3, takes currently months 
to reach the quasi-equilibrium state. As a result, it is not possible to obtain a full grid 
(spanned on BH mass, BH spin, and accretion rate) of disk solutions which might be 
applied for fitting to particular X-ray sources. Nevertheless, numerical modeling, either 
global or shearing-box, gives the insight into the real turbulent behavior of disk interior. 

The analytical, standard model of relativistic thin accretion disks (NOVIKOV & 
Thorne, 1973) has been in use for almost forty years. It successfully describes the 
disk component of the high/soft state spectra of AGN and microquasars. BH spins of a 
number of objects were estimated applying the NOVIKOV & Thorne (1973) model to 
their X-ray continuum spectra. Most recent global MHD simulations (KULKARNI ET 
AL., 2010) prove that, in the limit of thin disks, the emission and angular momentum 
profiles produced by MRI-generated viscosity are in good agreement with the predic- 
tions of the analytical, based on the ap prescription, model of NOVIKOV & Thorne 
(1973). 

However, as we have proven in the introduction to Chapter 4, some of the assump- 



174 



Conclusions 



tions behind the standard models of thin disks break down for accretion rates close to 
the Eddington value. Most importantly, the advective cooling triggers in and starts to 
affect disk observables. 

The model of slim accretion disks successfully generalizes standard thin disks to the 
regime of high luminosities. It accounts for the advective transport of heat and radial 
gradient of pressure which affects the profile of angular momentum. 

In this thesis we have presented the slim disk model in detail. We started from 
introducing the stationary, relativistic, height-averaged model (Chapter 4). We widely 
discussed the solutions, putting most attention to the locations of various disk inner 
edges. Basing on these solutions, a spectral model slimbb was constructed and de- 
scribed in Section 4.5.2. 

The initial model was generalized to the non- stationary case in Chapter 5. The 
equations describing time-evolution of a relativistic slim disk (Appendix B) have not 
yet been presented in the literature. We solved them using spectral methods and 
presented the limit cycle behavior of accretion disks around rotating BHs. We put 
most attention to the impact of BH spin on the disk observables. We showed that the 
only one significantly affected by BH spin was the maximal luminosity. We suggested 
a rapid, but very rough, method of estimating BH spin basing on the limit-cycle light 
curves. 

We developed a two-dimensional slim disk model including, in a self-consistent way, 
detailed treatment of the vertical radiative transfer in the disk interior (Chapter 6). 
We presented and discussed the vertical structure of disks and proved that the regular, 
one-dimensional models produce spectra which are very similar to the spectra obtained 
using the more precise, two-dimensional model. 

We also constructed a state-of-art model of self- irradiated slim disks (Chapter 7). 
The brief discussion included in the dissertation will be followed by a more detailed 
study to be published in a separate paper. 

In Chapter 8 we presented three possible applications of the slim disk models pre- 
sented in the thesis. First, we studied the spin evolution of BHs accreting at super- 
critical accretion rates. We proved that, under reasonable assumptions and for a rea- 
sonable set of disk parameters, the canonical spin value a* = 0.9978 (Thorne, 1974), 
may be exceeded. We also showed that for very high (m > 10) accretion rates, the evo- 
lution of BH spin is not influenced by the angular momentum and energy of photons 
captured by the BH. 

We applied the slim disk-based spectral model slimbb to high/soft spectra of LMC 
X-3. We suggested that the inconsistency in BH spin estimates for high luminosity data 
may be resolved if the viscosity parameter a decreases with luminosity. 

Finally, we also presented the way of normalizing global MHD simulations of accre- 
tion disks using the disk thickness profiles of slim disks. 

Slim disks are a powerful tool for modeling high/soft states of accretion disks at 
high luminosities. They can be further developed by including additional effects which 
may affect the real sources, and which are currently neglected or treated in a simplified 
way, e.g., outflows or modified viscosity prescriptions. As long as global, radiative MHD 
simulations of disks remain limited by the amount of computational power, slim disks 
will remain the reference model for luminous disks in AGN and X-ray binaries. 



Appendix A 
Kerr metric 



The specific angular momentum a of a BH witli mass M and angular momentum J is 
given by (G = c = 1), 

(A.1) 

Often a dimensionless parameter a* is introduced to express the magnitude of the spin, 

a* = — = — ^. (A.2) 

It ranges between a* = —1 (maximally counter-rotating BH), through a* = (Schwarzschild 
BH), to a* = 1 (maximally rotating BH). 

The space-time metric near rotating BHs was described by Kerr (1963). The 
general formula for the line element ds'^ is given in the Boyer-Lindquist coordinates 
[t,0,r,^] by, 

ds^ = -(i-'^)dt'-2'^asm'ed(f)dt+ (A.3) 
\ J P^ 

+ ^dr^ + p'de^ + (^^ + + sin^ o] sin' 



where 



A = - 2Mr + a^ (A.4) 

p' = r' + a' cos^ e. (A.5) 

The metric is stationary and symmetric with respect to the polar axis of ^ = and the 
equatorial plane 6 = 7r/2. The signature of the metric is ( — h ++)• 

For the purpose of studying accretion disks it is convenient to write the Kerr metric 
in the approximate form valid close to the equatorial plane (cos^ ^ ^ 0) and introduce 
the vertical cylindrical coordinate z = rcos6. The non-zero metric components are: 

^ _ r-2M ^ _ 2Ma „ _ A „_r-2 „_1 

9tt — — , 9t4> — — , 9<t><l) — :p2^ 9rr — 9zz — J-- 

tt _ A t<j) _ 2Ma 4,4) _ r-2M rr _ A_ zz _ i 

9 r^A' " rA ' " rA ' " ^.2i 9 -'-) 

where A = r'^ + r'^a? + 2Mra^. 
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A.l Killing vectors 

A vector field X is a Killing field if the Lie derivative with respect to X of the metric 
g vanishes: 

Cx9 = ^. (A.6) 

If the metric coefficients g^^, in some coordinate basis dx"" are independent of , then 
x^ = Sj^ is automatically a Killing vector, where Sj^ is the Kronecker delta (MiSNER 
ET AL., 1973). 

In the Kerr metric there are two obvious Killing vectors: 

connected with time and axial symmetries of the metric, respectively. Their contravari- 
ant components {[t,(j),r,6]) are, 

r/^ = (1,0,0,0), (A.9) 
e = (0,1,0,0), (A.IO) 

while the covariant take the form, 

Vt = V*9tt + V^9t<p = 9tu (A. 11) 

V<P = V^9t<t> + v'^9h = 9t^, (A. 12) 

^t = e9u + e9tr, = 9t^, (A.13) 
^4> = ^'9t4> + ^'^9H = 9h- (A- 14) 



Products of the Killing vectors yield. 





= rfvt = 


9tt, 




= ^% = 


9h 




= V% = 


9t<f>-, 




= i% = 


9t<f>' 



(A.15) 

(A.16) 
(A.17) 
(A.18) 



A. 2 Christoffel symbols 



The Christoffel symbols can be derived from the vanishing of the covariant derivative 
of the metric tensor gn; : 

= Vigik = - gmk^Ti - 9imX'M = - '^9m[k^i\i- (A. 19) 



A. 2 Christoffel symbols 
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By permuting the indices, and resumming, one can solve explicitly for the Christoffel 
symbols function of the metric tensor: 



^le - 2^"" + ^? ~ ^) - l^a'^^idrnk/ + 9me,k - 9ke,m), (A.20) 



Below we give the exact formulae for all the non-zero Christoffel symbols calculated for 
the general Kerr metric close to the equatorial plane (i.e., sin^ = cos^ 6 = 0). 
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Tl . TU ^ 'XiA (A.21) 



pep _ pep 
^ tr — ^ rt 



r 



pr 

pS _ p9 

r-e — Or 

pr pr 

J- — J- 



rr 



p</) _ p 

^ te — ^ 

pt pt 

r0 ~ 

r0 — 

■pi Y^t 



Ar2 

Ma 



Ar2 
2Ma2 cos 6* 

2Ma cose 



2\ 



Ma{3r'^ + a 

r<f> — (pr — 

r(r2 - 2Mr) - Ma^ 

rip ~ (pr ~ 

2Ma3 cos ^ 



rL = n. = (i + ^lcos^ 



2Ma2 

FL = TL = IH ^ cos a* 



MA 



MaA 

t(p — (pt — 



2Ma2 cos e 



2Ma(r2 + a2)cos^ 



P& _ py 

1 M - r 

r'^ = - + 

rr ' 



r A 
A 

r 

1 

r 

cos 6 



re - '-er - ^2 



cos 9 

a? cos ^ 

r2A 
A(Ma2 - r3) 

^4 



cos6'(Ar^ + 2Mr(r2 + a 



2\2\ 



j,6 



A. 3 Characteristic radii 
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The following characteristic radii, important for the accretion disk physics, may be 
distinguished (Bardeen et AL., 1972; Kato et AL., 2008): 
(i) The event horizon: 

The horizon is a boundary in spacetime beyond which events cannot affect an outside 
observer. Even light emitted from beyond the horizon can never cross it and leave the 
BH. The radius of the horizon is given by. 



M(i + v^r^). 



(A.22) 



(ii) Innermost stable circular orbit (marginally stable orbit): 
Defines the radius of the innermost stable circular geodesic orbit. It is located at. 



where 



nsco = M(3 + ^2 - V(3-Zi)(3 + Zi + 2Z; 

= 1 + (1 - aiy/-' 



[l + a. 



a/3 



a* 



1'3 



Sal + Zl 



(A.23) 

(A.24) 
(A.25) 



(iii) Marginally bound orbit 
Defines the innermost unstable circular orbit. Its radius is given by. 



2M (l - y + v^r^) . 



(A.26) 



(iii) Photon orbit 

At the photon orbit photons follow circular, unstable, trajectories around BH. The 
radius of such an orbit is, 



^ph 



2M 



1 + cos I - arccos(— a*^ 



(A.27) 



Fig. A.l presents locations of these characteristic radii versus dimensionless BH spin 
parameter a*. 
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Figure A.l: Characteristic orbits of the Kerr metric. The innermost stable circular 
orbit (risco)? the marginally bound orbit (rmb), the photon orbit (rph) and the BH 
horizon (rh) are presented. 



Appendix B 

Relativistic, non-stationary equations 
for slim disks 



B.l Metric and four- velocity components 

The Kerr metric coefficients close to the equatorial plane are given in Appendix A. 
The following form of the stress-energy tensor of the matter in the disk is assumed: 

= pu'u^ + pg'^ - t'^ + + u'q\ (B.l) 

where p is the rest mass density, is the viscous stress tensor and is the radiative 
energy flux. In the further part of the text we will also use the "vertically integrated" 
pressure P = 2Hp. 

The general form of the four-velocity is: 

u' = l {el,) + V^^^el^) + V^(^)e;^) + V^^^^eJ,)) , (B.2) 

where V^^^ are velocities as measured in the Local Non-Rotating Frame (LNRF) and 
cq) are LNRF basis vectors. Near the equatorial plane they take the following forms 
(BARDEEN ET AL., 1972): 

2Ma d 

AT/9. A IZ-J Q^' \^-^) 





^1/2 d 


Ht) = 








^(r) = 


r dr ' 




1 d 


6(6) = 


r do ' 




r d 


= 


AV2 dS' 
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According to Eq. B.2 the contravariant components of the four-velocity are: 

/11/2 

r 

, 2Ma , r 

'^1/2^1/2 ^ ^1/2' 

r 

We introduce the following velocity variables describing the flow: 

• the radial velocity of the fluid measured in the co-rotating frame (CRF) V (com- 
pare Eq. 4.4): 

, ^ I V 

7 V 1 - 

• the angular momentum per unit mass C: 

C = u^, (B.6) 

• the vertical velocity U in the form: 

cos 6. (B.7) 
By simple calculations we can express V^'^^ in terms of C: 

The Lorentz factor 7 at the equatorial plane is defined as following: 

7-2 = 1_ (00)2 ^gg) 

Taking into account Eqs. B.5 and B.8 we get: 

, / 1 £V-\ 1 £V ,„ , 

^'-[—. + —)-T^ + — (BIO) 

We assume that the Lorentz factor 7 at the disk surface is given by the same expression, 
as we expect the vertical velocities to be small compared to the radial and orbital 
motion. 

Now the contravariant components of the four-velocity take following forms: 

/11/2 

= (B.ll) 
V 



u'^ = — — + 7W 



= lu cose. 



B.2 Momentum conservation 
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The covariant four-velocity components are: 



Ut 
Ur 

ue 



-7 



AV2 

r V 



'yrU cos 6. 



[B.12) 



B.2 Momentum conservation 

Let us consider the r-component of the general equation: 

V,r'^ = ^ a" : 



1 ^„dp 1 f^2rpir\ 



-9 



[B.13) 
(B.14) 



dr Sr^ 

— h u — — 

ot or 

u^Tl^u' = ryu' + 2Tl^u'u^ + T^yu^ + Tl^u' + 2r;,«^«'' + Tl^ufut (B. 15) 

From among all components of the vertically integrated viscous stress tensor T**" we 
leave only T"' . Neglecting the terms proportional to cos^ 9 we get: 



u 



9t c?r 
What is equivalent to: 



p dr Sr^ dr 



(1 - V')grrU'— + — = - - + l^^^^A(^2y-) 

9t 1 — (9r r S 9r Sr^ (9r 

where A, Q, Q^, fi^, and R are defined as in Abramowicz ET AL. (1996): 



(B.16) 
(B.17) 



A 



MA 



1 - nm^ 



Q = 


u'P 2Mar 


A 


= 




^r3/2±aMV2 


Q = 


Q — U!, 


R = 


A 


r2AV2" 



7^3/2 



Finally: 



^3/2 



dt 7rAV2(i _ v^) 



V dV ^ A l-V^dp 1 - 1/2 a ^ ^ 



1 — V"^ dr r 



p dr EA dr 



ir'T"') 



(B.18) 
(B.19) 

(B.20) 
(B.21) 
(B.22) 

(B.23) 
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and 



dv_ _ yi - V^A 

'dt ~ 7^1/2 



V dV A l-V^dp d . 2^„, 



+ 



1 — dr r p dr SA dr 
The viscous stress tensor component T^^ is given by, 



[B.24) 



(B.25) 



where u is the kinematic viscosity coefficient and S is the surface density. The shear 
tensor components cXa/j are given in general by, 



1 1 



where B = m'^^ and h""^ = g"^ + u^u^ . 



(B.26) 



B.3 Angular momentum conservation 

Let us write the angular momentum conservation law in the form: 
where = S'^^y After vertical integration we get: 



(B.27) 



S u 



dt 



u 



dr 



- ^^{Tl) = 0, 



(B.28) 



where is vertically integrated (z, (f)) component of the viscous stress tensor. Let us 
follow Lasota (1994) and assume that the only nonvanishing component of Ti is: 



= 2i/S a^^ = z/L — . 



Finally: 



where: 



dt 



durh 1 d fui:A^/^A^/^-f^dn 



dr r dr 

2Mar r^A^/^C 



dr 



(B.29) 



[B.30) 



(B.31) 



A -lA^I'^ 

is the angular velocity with respect to the stationary observer. Introducing L and V 
we get: 



VA 



d£^ 

'dt ~ ~7V1 - y2AV2 97 ^ ^ 



dc AV2 d ^t^ sA^^/^Avy dn 

r^ dr 



(B.32) 



B.4 Disk stability criterion 
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B.4 Disk stability criterion 

Let us consider the vertical acceleration of the disk surface in the Zero Angular Momen- 
tum Observer (ZAMO) system of coordinates. The disk surface is defined in cylindrical 
coordinates as following: 

z = H{r,t). (B.33) 

Therefore, the infinitesimal shift in the vertical direction can be expressed in the fol- 
lowing way: 

dz = ^dr + ^dt. (B.34) 
or at 



Dividing both sides by dt we get: 

dz _ 
dt ~ 

Introducing U and V and substituting H = r cos Qh we obtain: 
d cos Qh U cos Qh 1 V d 



dr dt 



dt 



r'j \/l — V"^ dr 



r cos Qh] 



(B.35) 



(B.36) 



B.5 Vertical equilibrium 



Let us assume the following form of pressure (Abramowicz ET AL., 1997): 

ros^ f) 

P{r,e,t)=Po{r,t) 1 

L cos"^ Qh 

We will derive the vertical equilibrium equation taking 6 component of the 
equation: 

VkTt = 0. 

Simple manipulations give (after vertical integration): 

dP 



= Sn^VfcUe 



d9 ' 



1 dP 

^~de 



Since: 



cos 6 



[u^u^ - a^iutUt - 1)] 



dP 

'de 



2Pn 



-Oh 



COS Qh 



we finally obtain {for 6 = Qh and Pq = P): 
,duB 2P 



u 



dt 



E cos Qh 



, 2/ COS 6// ^.due _ 
+ [u^u^ - a [utUt - l)j — — u — = n, 



(B.37) 

general 
(B.38) 

(B.39) 
(B.40) 

(B.41) 
(B.42) 
(B.43) 

(B.44) 
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where: 



UtUt 



r- due 
Or 



u 



VI - Or 



{■yUr cos 6^) 



and 



u 



dt 



7- 



Ai/2 



dU U 



1- 



V dV Cr^ dC 



dt 7 V(l - V'^Y dt A dt 



Solving for ^ we get: 



Ai/2 



dU 



dt 7^1/2 cos 



T 



V dV Lr'^ dC 



(B.45) 
(B.46) 

(B.47) 



7f/ dcosQn 



cos 9 



H 



dt 



U dcosOn 



cos 9 



H 



dt 



(B.48) 



(B.49) 



where TZ is the right hand side of the Eq. B.44. 



B.6 Equation of continuity 



The general form of the continuity equation is: 

V,(pnO = 0, 

what can be written as: 

1 , . , 1 



V,(pu^) 



-9 



After the vertical integration it takes the following form: 



It*— - -- — (Era'') 
dt r dr 



dt • 



The time derivative ^ can be eliminated using the following relation: 



du^ 



^1/2 1 



dt rAV2 ^ 
Introducing the physical quantities we get: 

rAV2 



V dV Cr^ dC 

+ 



_{1-VY dt A dt 



'dt 



^AV2 



du' 19/ V Ai/2 
S— + -— I rS 



dt r dr \ y/Y^Ty^ r 



(B.50) 



(B.51) 



(B.52) 



(B.53) 



[B.54) 



B.7 Energy conservation 
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B.7 Energy conservation 

From the general form of the energy conservation equation: 

v,(rV) = 0, 



(B.55) 



we obtain, using the most general form of T*'^' (Eq. 2.4), in the non-relativistic fluid 
approximation , 



= p 



dt dr 



(B.56) 



where e is gas internal energy. Similarly like in LANDAU & LiFSHITZ (1959) (their Eq. 
49.4), after vertical integration, it could be put into the following form: 



ET 



dt 



dr 



vis /^rad 



g"^ - Q 



(B.57) 



Where S is entropy per unit mass, S is surface density of the disk, is the local 
viscous heat generation rate and Q'^^ is the radiative cooling rate Abramowicz ET AL. 
(1996): 



Q 



rad 



32aT^ 
3Sk 



Thermodynamical relations give (for p = p^ad + Pgas)'- 



dS 



^, IdT A dp 



T dr 
1 dT 



p dr 
Idp 



T dt ^ ^ ' pdt 



(B.58) 
(B.59) 

(B.60) 
(B.61) 



The appropriate sum of the derivatives of p = S/2i/ can be substituted using the 
continuity equation: 



f-ldp j.ldp du^ 1 d 2 r 



u -— — \- u -- 



p dt p dr dt dr 
Taking all together we obtain: 



— — rV + — . 



(B.62) 



Cv 

where 



dt dr 



( 



du* 1 9 , 2 r■^ 



Cv 



\ dt dr 

4-3/j P 
Eg - 1 ST' 

(4-3/3)(7g-l) 
12(l-/3)(7,-l)+/3' 



gvis_grad^ (B.63) 
(B.64) 

(B.65) 
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/3 = ^. 
P 



Ultimately: 

dT 
'dt 



1 

S 7^1/2 
VA 



g"^ - Q 

Cv 

dT 



rad 



7VI - F2v4V2 ar 



9m* 



^2 



[r u 



(B.66) 



(B.67) 
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Derivation of an arbitrary tetrad in 
the Kerr spacetime 

Our aim is to derive the tetrad of an observer moving along the photosphere that would 
depend only on the quantities which are calculated in accretion disk models, i.e., on 
the radial and azimuthal velocities of gas and the location of disk photosphere. 

The metric considered here is the Kerr geometry gik in the Boyer-Lindquist coor- 
dinates r, ^] (Appendix A). Similarly as in Carter's Les Houches lectures (Carter, 
1972), we will consider two fundamental planes; the symmetry plane Sq = [t,4>] and 
the meridional plane A^* = [r,6]. (Four)-vectors that belong to the plane Sq, will be 
denoted by the subscript 0, and vectors that belong to the plane A^*, will be denoted 
by the subscript *. For example, the two Killings vectors are ryg, Co- Note, that for any 
pair Xq, 17 one has, 

Xi,Y^'g,k = {XoY,)=0. (C.l) 



C.l Stationary and axially symmetric photosphere 
C.1.1 Photosphere 

Numerical solutions of slim accretion disks provide the location of the photosphere given 
by hpii{r) = rcosO. This may be put into rcos6 — /iph('") = F{r,6) = 0. The normal 
vector to the photosphere surface has the following [r, 6] components. 



dF OF 
dr ' 09 



N' 



dr 



where 



dg.(r) 
dr 



dF 

dr 



'dF 



(C.2) 



(C.3) 



is the derivative of the angle defining the location of the photosphere at a given ra- 
dial coordinate [cos^*(r) = h-pYi{f)/r\. Its non-zero components after normalization 
[{N,N,) = 1] are 
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grr f d6'* 
dr 



grr f 

gee V dr 



-1/2 



(C.4) 



1 -1/2 



There are two unique vectors S** confined in the [r, 9] plane which are orthogonal to 
A^* (and therefore are tangent to the surface). From [St^N^) = and (S^S^) = 1 one 
obtains the non-zero components of one of them: 



SI = (gr 



grr gee V 



-1/2 



(C.5) 



5 



. ^-l fde. 



grr gee V dr 



-1/2 



C.1.2 Tetrad 

The four-velocity Uq of an observer with azimuthal motion only is, 

4 = Ao iv' + ne) . (C.6) 

The normalization constant Aq comes from (moMq) = ~1 and equals 

Ao = [-gtt - ^g^iS^ - 2a;)]-'/' . (C.7) 

It is useful to construct a spacelike vector (kq) confined in the [t, 0] plane, that is 
perpendicular to Uq. From (kk) = 1 and (kmq) = we have. 



1 -fi/)(l -w/)] 



1/2 • 



(C.8) 



where / = —Urp/ut is the specific angular momentum and u = —gt(t,/g4>4i is the frequency 
of frame dragging. The set of vectors [mq, A^*, Kq, S]] already forms the desired tetrad 
valid for the pure rotation {u^ = 0) case. 

Similarly, the four-velocity u of gas moving along the photosphere with non-zero 
radial velocity may be decomposed into 



u 



A{^^ + QC + vS:). (C.9) 
The normalization condition (uu) = —1 gives, 

A = [-gu - ^gu{^ - 2uj) - v^] , (C.IO) 
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where v is related to the radial component of the gas four-velocity u*" by: 



' i ur/S:f i-gu - ng^^n - 2a;)) 
l + {ur/S:f 

Taking into account that 5*^ and are the only non-zero components of vector 5"*, it 
is obvious that is orthogonal to Kq. 

The vectors we have just calculated {u, kq) are both orthogonal to A^* since (A^*5'*) = 
0. To complete the tetrad we need one more spacelike vector [S) that is orthogonal to 
these three. Let us decompose it into, 

= au' + PKl + -fNi + 6Si. (C.12) 

The orthogonality conditions {kqS) = and (N^S) = give immediately 7 = /3 = 0. 
The only non-trivial condition is (uS) = 0. Together with (SS) = 1 it leads to: 

= (1 + iV)-^/2(ii;n* + Si). (C.13) 
The vectors m*, A^*, Kq, form an orthonormal tetrad in the Kerr spacetime: 

e^(^) = [u\ Nl, 4, S% (C.14) 

This tetrad is known directly from the slim disk solutions, as it depends on the calcu- 
lated quantities (u'', Q, I and 0^,{r)) only. Any spacetime vector X*, could be uniquely 
decomposed into this tetrad with, X(^) = Xje*^^^. 



C.2 General case 

In this section we will assume nothing about the four-velocity of matter and the loca- 
tion of photosphere. Both may be non-stationary and non-axially symmetric. Following 
the same framework as in the previous subsection, we will describe how to obtain the 
tetrad of an observer instantaneously located at the photosphere that depends only on 
the quantities calculated by accretion disk models. 



C.2.1 Photosphere 



In the most general case of a non- stationary and non-axially symmetric photosphere, 
the location of the photosphere may be described by the following condition. 



F(t,0,r,^)=O. 
The vector A^ normal to the photosphere has the components. 



OF OF OF OF 
dt ' d(h' dr ' 89 



(C.15) 



(C.16) 



which may be calculated from slim disk solutions. 
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Derivation of an arbitrary tetrad in the Kerr spacetime 



Let us project N into the instantaneous 3-space of the comoving observer (C.22) and 
normahze to a unit vector after the projection, 

N' = ^Jt , N' = N^h\. (C.17) 

|(iVAr)|i/2' 

In terms of the tetrad C.24, such constructed vector A^* has the following decomposition, 

N^ = N[a (n^) + 1 {Nl) + 7 K) + 5 {8% (C.18) 
The components iV, a, 7, 5 are known. 

C.2.2 Tetrad 

Similarly as in Section C.1.2, we may always uniquely decompose u*, a general timelike 
unit vector, into 

u' = A{v^, + vS:), (C.19) 

where Uq is a timelike unit vector, and SI is a spacelike unit vector. Formula C.19 
uniquely defines the two vectors Uq^SI and the two scalars A, v. The vectors and 
scalars 

{Avy,,si}, (C.20) 

can be calculated from known quantities given by slim disk model solutions. 

The four-velocity (Eq. C.19) defines also the instantaneous 3-space of the comoving 

observer with the metric 7jfc and the projection tensor /i^*, 

lik = 9ik~UiUk, (C.21) 
h\ = 5\-u'uk (C.22) 

We define two unit vectors Hq and A'^^ by the unique condition, 

(kqUo) = 0, {S,N,) = 0. (C.23) 

As before, the four vectors 

e\^) = K, 4, A^:, S% (C.24) 

which form an orthonormal tetrad of an observer with the four- velocity Uq can calculated 
from the solutions of the slim disk equations. 

Let us now write decompositions of the four vectors: the first two we have derived, 
the next two guessed (but the guess should be obvious): 

= A[l{ul) +0{N:) +0{kI) +V{S% (C.25) 

= N[a(ul) + im +6{S:)], (C.26) 

= k [0{ul) +b{Nt} +1(4) +0(5:)], (C.27) 

= S[A{uI,) + B{N:) + C{kI) + 1{S:)]. (C.28) 
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The four unknown components, b,A,B,C one calculates from the four non-trivial or- 
thogonality conditions {{un) = by construction, cf. Eqs. C.25 and C.27), 

(uS) = 0, (NS) = 0, (Sk) = 0, (Nk) = 0, (C.29) 

and the two unknown factors k, S, from the two normalization conditions 

(kk) = -1, {SS) = -1. (C.30) 

The conditions C.29, C.30 are given by linear equations. 

Equations C.25-C.30 define the tetrad e of an observer conioving with matter, 
and instantaneously located at the photosphere: 

e = [u\ N\ k\ S']. (C.31) 

Both the matter and the photosphere move in a general manner. The zenithal direction 
in the local observer's sky is given by A^*. 
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Appendix D 

Integration over the world-tube of the 
photosphere 



For stationary and axially symmetric models we have so far (Appendix C): 

= Nl = unit vector orthogonal to the photosphere. It is in the [r, 6] plane, 

SI = unit vector orthogonal to A^* that lives in the [r, 6] plane, 

= four- velocity of matter. It lives in the [t, 0, r, 6] space-time, 

unit vector orthogonal to f/' that lives in the [t, 0] plane, 

S"* = unit vector orthogonal to W, N''' and /t*. It lives in the [t, 0, r, 6'] space-time, 

e = [u*, A^*, K*, 5*] = the tetrad comoving with an observer located in the 
photosphere. 

The integration of a vector over the 3-D hypersurface orthogonal to A^* (i.e. the 
3-D world-tube of the photosphere) may be symbolically written as 



where dS is the "volume element" in "H. 

Obviously, the hypersurface "H is spanned by the three vectors [m*, k*, S'^]^. Each of them 
is a linear combination of [77*, Sl]]sf, and each of the three vectors from [77*, SI]n is 
orthogonal to A^*. 

Therefore, one may say that the hypersurface "H is spanned by [ri\C,\ SI]n. It will be 
convenient to write 



where dR is the line element along the vector SI, i.e. along the photosphere in the [r, 6] 
plane, with 6 = 6^,{r) defining the location of the photosphere, and dA is the surface 
element on the [t, 0] plane. 




(D.l) 




(D.2) 
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Integration over the world-tube of the photosphere 



In order to calculate dA, imagine an infinitesimal parallelogram with sides that are 
located along the t = const and = const lines. The proper lengths of the sides are 
du = \gtt\^^'^dt and dv = Ig^^l^^'^dcp respectively, and therefore dA, which is just the 
area of the parallelogram, is given by 

dA = du dv sin a = dt d0 \gtt\^^'^\g(i,4,\^^'^ sin a, (D.3) 

where a is the angle between the two sides. Obviously, the cosine of this angle is given 
by the scalar product of the two unit vectors rii and Xi pointing in the [t, </>] plane into 
t and (j) directions respectively. These vectors are given by (note that rij = ZAMO), 

n. = m m (D4) 

Because (Vjt) = 5*^ and (Vi0) = 5*^^, one may write, 

^ |^ti|l/2' ^« ^ |^00|l/2- ^^-^^ 

Therefore, 
and 

Inserting this into the formula for dA we get dA = dtd(f){g^^ — gttg(j>(i>Y^'^ ■ The final 
formula for dS* is. 



dS = dt d0 dr {g^^ - gu g^^) ^'"^ J grr + gee (^^^ • (D-8) 
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